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Chapter  1 


Introduction 


Turbulent  fluid  motions  are  typically  characterized  by  several  features  including  randomness 
in  both  space  and  time,  vorticity,  an  energy  cascade  from  large  to  small  scales  where  energy 
dissipation  occurs,  and  a  large  increase  in  diffusion  of  properties  (he.,  temperature,  salinity) 
compared  with  molecular  diffusion  [McDougal  et  ah,  1988],  These  features  of  turbulent  flows 
are  usually  caused  by  some  sort  of  flow  instability.  In  homogeneous  flows,  instabilities  are 
usually  related  to  the  Reynolds  number,  which  can  be  thought  of  as  the  ratio  of  inertial  to 
viscous  forces.  As  Reynolds  number  increases,  inertial  forces  overcome  viscous  dissipation, 
and  instabilities  grow  until  they  overtake  the  flow.  In  a  density  stratified  flow,  a  gravity 
force  is  present  which  acts  as  a  stabilizing  force,  giving  rise  to  a  buoyancy  force  that  must 
be  overcome  as  well  as  viscous  forces  for  the  fluid  to  become  turbulent.  This  suggests  that 
the  Fronde  number,  which  can  be  considered  a  ratio  of  inertial  to  gravity  forces,  will  also  be 
important  in  determining  if  a  flow  becomes  turbulent.  Moreover,  the  addition  of  a  buoyancy 
force  helps  create  a  situation  where  the  transition  to  turbulence  is  marked  by  intermittent 
turbulent  patches  in  the  flow,  rather  than  a  smooth  transition  throughout  the  flow  as  in  a 
homogeneous  flow.  The  characterization  of  this  intermittcncy  of  turbulence  within  a  density 
stratified  flow  is  an  area  of  active  research. 

When  one  considers  the  fact  that  the  atmosphere  and  ocean  are  density  stratified  fluids, 
it  can  be  said  that  the  vast  majority  of  the  flows  on  Earth  take  place  in  stratified  fluids. 
They  occur  in  the  ocean  below  the  mixed  layer,  in  the  stratosphere,  and  in  the  nocturnal 
atmospheric  boundary  layer  (at  night  the  sun  does  not  provide  energy  to  mix  the  atmosphere 
near  the  Earth).  Turbulence  often  occurs  in  these  flow  regimes  due  to  free  shear  instabilities 
(e.g.,  Kelvin-Helrnholtz),  internal  wave  breakdown,  and  wakes  of  structures  such  as  islands, 
mountains,  and  submarines.  Turbulence  in  these  areas  has  impact  ranging  from  weather 
prediction  to  pollution  dispersion.  It  is  this  area  of  geophysical  turbulence  that  will  be  the 
focus  of  this  dissertation. 

An  important  effect  of  stratification  is  that  gravity  allows  internal  gravity  waves  to 
form.  Internal  waves  have  the  ability  to  propagate  energy  throughout  the  flow  [e.g.,  Riley 
and  Lelong,  2000,  Slinn  and  Riley,  1996,  Winters  et  al.,  1995].  Turbulence  caused  by  the 
breakdown  of  internal  gravity  waves  can  affect  the  mixing  of  heat  and  elements  within 
the  fluid  [Lombard  and  Riley,  1996a, b,  Slinn  and  Riley,  1996,  1998].  While  important  in 
density  stratified  flows,  breakdown  of  internal  waves  will  not  be  examined  in  this  study. 
Rather,  the  focus  of  this  study  will  be  turbulence  formed  by  free  shear  instabilities  by 
flows  dominated  by  vortical  modes.  The  remainder  of  this  chapter  contains  description 
of  density  stratification  that  occurs  in  nature,  followed  by  a  literature  review  of  density 
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stratified  flows.  Chapter  2  contains  theoretical  considerations  including  a  derivation  of  the 
equations  of  motion  and  a  description  of  the  Boussinesq  approximation.  Chapter  3  contains 
results  of  studies  that  investigate  vertical  shear  and  dissipation  rate,  buoyancy  Reynolds 
number,  and  parameterization  of  turbulence  using  Taylor-Green  simulations  and  uniform 
density  stratification.  Chapter  4  contains  results  of  investigations  that  study  the  effect 
of  non-uniform  density  gradient  on  simulated  wake  flow,  including  a  method  of  calculating 
potential  energy  in  non-uniform  stratification.  Chapter  5  contains  conclusions  and  proposed 
future  work. 

Notation  Several  equations  are  introduced  in  this  document,  and  often  it  is  difficult  to 
distinguish  between  dimensional  and  nondimensional  quantities.  As  such,  a  convention  is 
adopted  in  this  document  where  (:)  is  used  to  denote  dimensional  quantities,  while  unmarked 
quantities  (i.e.,  no  tilde)  will  denote  nondimensional  quantities. 


1.1  Density  Stratification  in  Fluids 

In  the  atmosphere  and  ocean,  distinct  layers  form  that  are  characterized  by  the  rate  of 
temperature  (and  hence  density)  change  with  height.  The  atmosphere  is  typically  separated 
into  four  layers;  starting  from  the  Earth’s  surface  and  increasing  in  height  these  layers  are 
the  troposphere,  stratosphere,  mesosphere,  and  thermosphere  [e.g.,  Brasseur  and  Solomon, 
1984,  Labitzke  and  van  Loon,  1999,  Lutgens  and  Tarbuck,  1995].  The  ocean  is  typically 
divided  into  three  regions;  starting  from  the  ocean  surface  and  increasing  in  depth  these 
layers  are  the  mixed  layer,  the  thermocline,  and  deep  water  [Colling,  2002,  Gill,  1982]. 
This  section  contains  a  brief  description  of  each  layer  and  its  significance  to  geophysical 
turbulence,  and  the  reader  is  referred  to  the  aforementioned  references  for  further  discussion 
on  the  atmosphere  and  ocean. 


1.1.1  Atmosphere 


The  troposphere  is  the  lowest  layer  in  the  atmosphere.  It  varies  in  height,  spanning  from  the 
Earth’s  surface  to  approximately  18km  over  the  tropics,  while  spanning  to  approximately 
10km  over  the  Earth’s  poles.  The  troposphere  contains  approximately  80%  of  the  total  air 
mass  of  the  Earth,  and  is  where  all  weather  phenomenon  takes  place.  During  the  day  solar 
heating  of  the  Earth’s  surface  give  rise  to  convection  currents  in  the  troposphere,  causing 
the  troposphere  to  be  well  mixed.  At  night  solar  heating  ceases,  and  radiative  cooling  causes 
a  stable  density  stratification  layer  to  form  in  the  lower  3-5kin  of  the  troposphere.  Hence, 
this  layer  is  referred  to  as  the  nocturnal  boundary  layer. 

The  adiabatic  change  in  temperature  with  height  is  called  the  lapse  rate,  G,  and  is 
defined  as  [Gill,  1982,  p.50]: 


G  =  gotT /  cp , 


(1.1) 

(1.2) 


where  g  is  the  gravitational  acceleration,  T  is  the  temperature,  and  cp  is  the  specific  heat, 
5  is  the  thermal  expansion  coefficient  at  constant  pressure  /5,  and  humidity  (or  salinity  in 
the  ocean)  5: 


_  _  1  dp 

POT 


p,S 


(1.3) 
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Since  there  is  a  negative  sign  in  front  of  G  in  (1  1),  and  altitude  is  measured  from  the  ground 
upwards,  G  is  a  measure  of  the  decrease  in  temperature  with  height.  G  can  (and  does) 
become  negative  (signifying  an  increase  in  temperature  with  altitude),  particularly  in  the 
stratosphere  and  thermosphere.  G  of  the  troposphere  is  taken  on  average  to  be  6.5°C/km, 
but  can  vary  locally  and  depends  on  humidity  content  .  For  example,  the  nocturnal  boundary 
layer  is  stably  stratified,  and  has  a  negative  G.  Also,  it  is  possible  for  cooler  air  to  be  trapped 
near  the  Earth’s  surface,  resulting  in  a  locally  stable  stratification  layer  with  positive  G. 
An  area  of  cooler  air  trapped  near  the  surface  is  termed  a  thermal  inversion,  and  can  cause 
pollutants  and  smog  to  be  trapped  near  the  Earth’s  surface  instead  of  being  convected  away 
from  the  surface.  Thermal  inversions  commonly  occur  when  cool,  moist  air  from  the  ocean 
blows  over  land  and  when  a  warm  front  moves  into  a  region,  trapping  low  temperature 
underneath.  Westerly  winds  off  the  Pacific  Ocean  make  Los  Angeles  an  ideal  location  for 
temperature  inversions  to  form,  causing  the  city’s  famous  smog. 

The  upper  bound  of  the  troposphere  is  the  tropopause.  The  tropopause  is  marked  by  a 
sharp  change  in  G,  and  signifies  the  boundary  between  the  troposphere  and  stratosphere. 
The  location  of  the  tropopause  decreases  in  altitude  from  approximately  18km  in  the  tropics 
to  approximately  10km  near  the  poles. 

The  stratosphere  spans  from  the  tropopause  to  approximately  50km  above  the  Earth’s 
surface,  and  with  the  troposphere  contains  99%  of  the  Earth’s  air  mass.  The  lower  2  10km 
of  the  stratosphere  is  near  isothermal  (i.e.,  G  is  near  0°C/km),  while  above  20km  G  becomes 
negative  (increase  in  temperature  with  height).  The  negative  G  is  believed  to  be  caused  by 
the  absorption  of  ultraviolet  radiation  from  the  sun  by  ozone.  Since  density  decreases  with 
temperature  for  most  (if  not  all)  gases,  the  increase  in  temperature  with  height  causes  a  very 
stable  density  stratification  (hence  the  name  stratosphere).  The  stable  stratification  acts  as 
a  barrier  to  vertical  currents  from  the  troposphere  and  inhibits  vertical  motion  within  the 
stratosphere.  Such  inhibition  of  vertical  motion  will  result  in  motions  to  preferentially  grow 
in  the  horizontal,  which  is  of  interest  regarding  mixing  of  elements  and  chemical  reactions 
that  take  place  in  the  stratosphere. 

The  area  where  G  becomes  positive  again  is  called  the  stratopause,  and  marks  the 
separation  between  the  stratosphere  and  the  mesosphere. 

The  mesosphere  and  thermosphere,  while  playing  vital  roles  in  heat  absorption  from 
the  sun,  are  uninteresting  from  a  fluid  mechanics  perspective.  The  mesosphere  ranges 
from  approximately  50km  to  85km.  while  the  thermosphere  ranges  from  85km  to  500km. 
The  lapse  rate  in  the  mesosphere  is  approximately  2°C/km,  while  in  the  thermosphere 
temperature  increases  with  height  as  solar  photons  are  absorbed.  It  is  interesting  to  note 
that  the  northern  lights  occur  in  the  thermosphere,  and  that  the  international  space  station 
has  a  stable  orbit  in  the  upper  thermosphere 

1.1.2  Ocean 

The  mixed  layer  of  the  ocean  varies  in  depth,  spanning  from  the  ocean  surface  to  a  depth 
of  approximately  10m  near  the  poles,  200m  in  the  mid-latitudes,  and  50m-100m  in  the 
tropics.  This  layer  is  called  the  mixed  layer  because  wind  and  waves  cause  it  to  be  well 
mixed,  resulting  in  temperature  and  salinity  profiles  that  are  close  to  uniform.  Nearly  all 
sunlight  is  absorbed  in  the  mixed  layer,  causing  temperatures  in  excess  of  30°C. 

The  lapse  rate  G  of  seawater  is  much  smaller  than  that  of  the  atmosphere,  since  air  is 
much  more  compressible  that  seawater.  G  has  a  typical  value  of  0.125°C/km,  and  ranges 
between  0.1  and  0.2°C/km  [Fofonoff  and  Millard,  1983,  p.38].  Note  that  G  is  positive  for 
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the  ocean,  but  marks  decrease  in  temperature  with  depth,  since  depth  is  measured  from 
sea  level  downward  from  the  sea  surface  (as  opposed  to  the  atmosphere  where  height  is 
measured  upward  from  the  surface). 

The  thermocline  ranges  from  the  bottom  of  the  mixed  layer  to  approximately  1000m. 
The  region  is  called  the  thermocline  because  temperature  drops  rapidly,  as  much  as  40°C/km 
near  the  top  of  the  layer,  cooling  the  thermocline  to  near  10°C  at  its  lower  boundary.  This 
drop  in  temperature  causes  an  increase  in  density  with  depth,  hence  the  region  is  sometimes 
referred  to  as  the  pycnocline.  Since  density  increases  with  depth,  the  thermocline  is  stably 
stratified.  Since  the  thermocline  is  a  link  between  the  surface  and  the  deep  ocean,  the 
dynamics  of  the  thermocline  are  important  and  play  a  role  in  such  topics  as  pollution 
dispersion,  motion  of  small  food  sources  (e.g.,  plankton),  and  heat  dissipation  between  the 
surface  and  deep  ocean. 

The  deep  ocean  extends  from  the  bottom  of  the  thermocline  to  the  ocean  floor.  In 
the  deep  ocean  temperature  and  salinity  (and  hence  density)  are  relatively  uniform.  In¬ 
terestingly,  local  patches  of  high  density  water  are  formed  at  the  surface  in  cold  regions 
including  the  North  Atlantic  off  the  coast  of  Greenland  and  in  the  Antarctic  near  the  Ross 
and  Weddell  Seas  where  surface  water  freezes.  When  ocean  water  freezes  it  is  fresh  water 
that  converts  to  solid  mass,  leaving  a  higher  salinty,  higher  density  fluid.  (This  process 
is  called  brine  rejection  [Colling,  2002,  p.213]).  These  cold,  high  density  sources  of  water 
sink  from  the  surface  and  through  the  thermocline,  forming  a  convective  current.  This 
convective  current,  along  with  wind  driven  currents,  is  referred  to  as  the  “conveyor  belt,” 
and  helps  bring  water  from  the  deep  ocean  to  the  surface  in  the  mid-latitude  Pacific.  The 
rate  of  travel  along  the  conveyor  belt  is  on  the  order  of  1000  years. 


1.2  Review  of  Previous  Work 


1.2.1  Wake  turbulence  in  stratified  fluids 


Wakes  are  generated  when  there  is  relative  motion  between  a  body  and  adjacent  fluid. 
Wakes  in  stratified  fluid  can  be  found  in  many  settings,  including  such  as  those  generated 
by  airplanes,  and  submarines,  mountains,  and  buildings.  When  analyzing  wakes  in  stratified 
fluids,  the  Froude  (F)  and  Reynolds  (Re)  numbers  are  typically  defined  in  terms  of  the  object 
size  and  velocity.  Many  experiments  involve  towing  a  sphere  in  a  tank  containing  a  density 
stratified  fluid,  leading  to  common  F  and  Re  definitions: 


F 


U  UR 

— ,  Re=— , 
NR  v 


where  U  is  sphere  velocity,  R  is  sphere  radius,  and  N  =  g/podp/dz  is  buoyancy  frequency. 
Here  g  is  the  gravitational  acceleration  constant,  po  is  a  reference  density,  and  p  is  the 
background  density. 

Scaling  arguments  from  Spedding  et  al.  [199Ga]  show  that  F  >  3  for  initial  active  turbu¬ 
lence  (turbulence  with  the  ability  to  overturn)  to  occur.  Chomaz  et  al.  [1993]  demonstrate 
that  the  wake  behaves  initially  as  a  homogeneous  fluid  for  F  >  4.5.  In  the  late  wake,  many 
researchers  have  found  these  initial  vortices  to  increase  in  horizontal  length  and  decrease  in 
vertical  length,  forming  “pancake  eddies”  [Bonnier  et  al.,  1998,  Flor  et  al.,  1995,  Spedding 
et  al.,  199Ga,b].  In  addition,  Spedding  et  al.  [199Ga,b]  have  shown  the  wake  width  profile 
to  be  the  same  as  an  unstratified  flow,  but  with  peak  velocity  an  order  of  magnitude  larger 
than  unstratified  flow. 


1.2 .  REVIEW  OF  PREVIOUS  WORK 


It  has  been  noted  that  horizontal  planes  in  the  wake  are  similar  to  monopole  or  dipole 
formation  [Riley  and  Lelong,  2000].  While  many  studies  on  the  evolution  of  dipoles  have 
been  carried  out  [Garten  et  ah,  1998,  Fraud  and  Fincham,  2003,  Spedding,  2002],  an  in¬ 
teresting  study  is  that  of  Billant  and  Cliomaz  [2000b].  They  show  that  vertical  columnar 
dipoles  undergo  a  “zig-zag”  instability,  causing  the  dipoles  to  be  divided  into  separate  pan¬ 
cake  vortex  layers.  They  note  this  instability  to  occur  between  0.13  <  F^o  <  0.21,  with 
Fho  =  Uo/(NR)  based  on  the  initial  dipole  traveling  velocity  and  radius. 

1.2.2  Mixing  Efficiency 

Mixing  is  a  small  scale  process  affecting  the  thermodynamic  makeup  of  a  fluid.  It  is  ir¬ 
reversible,  as  the  fluid  can  not  be  returned  to  its  original,  pre-  mixed  state.  Mixing  is 
important  in  geophysical  flows  as  it  relates  directly  to  the  dynamics  of  heat,  chemicals,  and 
pollutants  in  the  atmosphere  and  ocean.  Mixing  is  often  quantified  by  an  “efficiency”,  or  a 
relation  between  the  rate  of  conversion  of  available  potential  energy  to  background  potential 
energy  (each  defined  below)  to  the  rate  of  lost  to  viscous  dissipation.  There  are  differing 
definitions  of  mixing  efficiency  in  the  literature  which  will  be  examined  in  the  following 
paragraphs.  Prior  to  defining  mixing  efficiency,  the  concepts  of  available  and  unavailable 
potent  ial  energy  are  discussed. 

Discussion  of  potential  energy  in  geophysical  settings  usually  involves  the  concepts  of 
available  and  background  potential  energy,  first  suggested  by  Lorenz  1955].  He  noted  that 
in  order  to  convert  the  total  potential  energy  in  the  Earth’s  atmosphere  to  kinetic  energy 
the  temperature  needed  to  reach  absolute  zero,  and  all  mass  would  be  located  at  sea  level 
(z  =  0);  conditions  that  cannot  readily  occur.  (It  is  estimated  that  potential  energy  makes 
up  25%  of  the  total  energy  (internal  +  potential  +  kinetic)  in  the  Earth’s  atmosphere, 
while  only  2%  is  kinetic  energy  [Gill,  1982,  pg.81]).  Instead,  the  potential  energy  Ev  that 
is  available  for  conversion  to  kinetic  energy  is  said  to  be  the  result  of  any  deviation  from  a 
background  (or  rest)  potential  energy,  Eb.  Eb  is  a  state  that  would  exist  if  the  fluid  were 
adiabatically  redistributed  (i.e.,  no  heat  transfer)  to  a  minimum  energy  state.  The  available 
potential  energy  is  the  total  potential  energy,  V ,  minus  the  background  potential  energy: 

EP  =  V  -  Eb. 

Here  V  is  defined  as  p(z)gz,  and  where  p(z)  is  density  as  a  function  of  vertical  position  z 
and  g  is  the  gravitational  acceleration.  Due  to  fluid  motions,  at  a  given  instant  p(z)  may 
not  be  in  its  lowest  energy  state.  Eb  is  the  minimal  potential  energy  attainable  through 
the  adiabatic  redistribution  of  p  [Thorpe,  1977,  Winters  et  ah,  1995]  (further  discussion  on 
the  computation  of  Eb  is  found  in  §4.2.2),  and  any  change  to  Eb  is  deemed  mixing.  Similar 
explanations  for  Eb  and  Ea  can  be  found  in  Staquet  [2000]  and  Peltier  and  Caulfield  [2003]. 

There  are  two  definitions  of  mixing  efficiency  in  the  literature.  One  definition  of  mixing 
efficiency  is  the  ratio  of  energy  lost  to  background  potential  energy  to  the  rate  of  kinetic 
energy  dissipation  lost  to  internal  energy  [e.g  Winters  et  al.,  1995]: 

r«  =  |,  (i.4) 

where  \  is  the  irreversible  rate  of  potential  energy  dissipation  to  background  or  unavailable 
potential  energy.  A  second  definition  of  mixing  efficiency  is  the  proportion  of  kinetic  energy 
lost  by  the  fluid  that  leads  to  mixing,  leading  to  the  relation  [e.g.,  Peltier  and  Caulfield, 
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2003]: 


r6  =  ^. 

X  +  e 


(1.5) 


Both  ra  and  T are  easily  computed  in  numerical  simulations.  However,  in  field  ex¬ 
periments,  x  and  c  are  difficult  to  make  at  the  same  time  due  to  time  requirements  for 
measuring  each  quantity  [Gargett  and  Mourn,  1995].  It  would  he  beneficial  to  measure  one 
quantity  (usually  e)  and  relate  it  to  the  other.  This  is  done  through  the  flux  Richardson 
number  Rf,  defined  as  either  the  fraction  of  shear  turbulent  kinetic  energy  spent  increasing 
potential  energy  [Smyth  et  al.,  1990]  or  the  ratio  of  loss  of  kinetic  energy  to  buoyancy  flux 
to  that  produces  by  shear  [Gargett  and  Mourn,  1995].  Both  definitions  lead  to  the  relation: 


Rr  = 


B 


B+I 


(1.6) 


where  B  =  —  (g/pQ)p'w;  is  the  turbulent  buoyancy  flux,  and  primes  indicate  fluctuating 
component  (with  w'  vertical  fluctuating  velocity).  In  its  current  form,  (1.6)  does  not  relate 
to  mixing,  just  the  fraction  of  the  loss  of  kinetic  energy  put  into  potential  energy.  To  relate 
this  to  mixing,  density  is  assumed  to  be  linear  in  temperature,  and  the  buoyancy  flux  can 
be  written  in  terms  of 


0T\20p 
Oz  )  dz 


This  allows  Rf  to  be  related  to  mixing  by  means  of  a  dissipation  flux  coefficient: 


Th  = 


Rf 


1  —  Rf 

and  mass  diffusivity  parameter  [Osborn,  1980]: 


Kp  —  I T  — 

P  N2 


(1.7) 


(1.8) 


In  oceanic  applications,  is  often  considered  a  constant  with  value  0.2  as  suggested  by 
Osborn  [1980].  This  value  is  then  used  to  compute  \  based  on  measures  of?  [Smyth  et  ah, 
1996].  However,  measurements  of  heat  flux  and  ?  have  shown  to  vary  between  0.1  to  0.4 
in  the  open  ocean  [Mourn,  1990]  to  0.7  in  turbulent  tidal  front  [Gargett  and  Mourn,  1995]. 


1.2.3  Fossil  Turbulence 

The  topic  of  fossil  turbulence  is  championed  by  Gibson  [1980],  although  he  states  the  con¬ 
cept  has  been  developed  by  others  as  early  as  1969.  Fossil  turbulence  is  defined  as  “the 
remnant  fluctuation  in  any  hydrophysical  field  produced  by  active  turbulence  which  persists 
after  the  fluid  is  no  longer  actively  turbulent  (overturning)  on  the  scale  of  the  fluctuation” 
[Gibson,  1980,  1987].  Gibson’s  hypothesis  is  that  the  vast  majority  of  the  ocean  is  fossilized 
turbulence  with  very  few  patches  of  active  turbulence. 

The  concept  of  fossilization  is  as  follows:  an  “actively  turbulent”  patch  is  formed  by  some 
instability.  In  an  actively  turbulent  patch,  inertial  forces  are  greater  than  buoyancy  forces, 
and  vertical  overturning  will  occur.  The  onset  of  fossilization  occurs  when  the  inertial  forces 
of  the  largest  scales  become  equal  to  the  buoyancy  forces.  In  this  state,  the  vertical  length 
scale  of  the  patch,  Lp,  becomes  equal  to  the  overturning  scale  L0  (see  (2.53)).  These  large 
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scales  are  now  considered  fossilized  since  they  can  not  overturn.  As  energy  is  dissipated  (via 
?),  the  inertial  force  of  the  large  scales  diminishes.  Then,  without  energy  input,  buoyancy, 
inertial,  and  viscous  forces  will  all  be  equal.  At  this  state,  assuming  isotropy,  a  transition 
kinetic  energy  dissipate  rate  is  defined  as  7tr  =  257'N2. 

Once  fossilized,  the  patch  is  now  in  a  state  of  microstructure  events.  The  fossilized 
patch  then  exhibits  “secondary  turbulence  events,”  presumably  due  to  instabilities  within 
the  original  turbulent  patch.  These  secondary  turbulence  events  entrap  additional  fluid  into 
the  fossilized  patch.  The  size  of  the  patch  will  expand,  and  then  fossilize  by  the  pathway 
described  above.  These  secondary  events  are  used  to  explain  large  regions  of  microstructure 
activity  measured  in  the  ocean. 

The  concept  of  fossilization  is  not  widely  agreed  upon.  Gregg  [1987]  cites  calculations  of 
time  constants  that  suggest  the  maximum  possible  age  of  the  microstructure  is  well  below 
that  required  for  fossilization.  Also,  Gregg  points  to  assumptions  made  by  Gibson,  namely 
that  ?  (and  mixing)  are  maximum  at  turbulence  collapse,  while  measurements  suggest 
significant  mixing  after  turbulent  collapse. 

1.2.4  Mixing  in  Ocean  Boundaries 

Vertical  mixing  in  the  ocean  is  important  to  the  transport  of  heat  and  salt  as  well  as 
supplying  nutrients  from  bottom  water  to  the  surface  for  support  of  biological  life.  Vertical 
mixing  is  often  related  to  an  “eddy  diffusivity,”  which  is  the  value  used  to  describe  the  mixing 
caused  by  eddy  motion.  In  the  open  ocean,  measured  values  of  vertical  eddy  diffusivity  are 
approximately  10-20%  of  what  is  required  for  vertical  mixing  to  occur,  suggesting  that  very 
little  vertical  mixing  occurs  in  the  open  ocean  [Slinn  and  Riley,  1996].  Instead,  this  results 
suggests  that  the  majority  of  mixing  in  the  ocean  occurs  at  the  ocean  boundaries  (e.g.. 
near  continental  slopes,  islands,  and  other  topological  features).  One  method  of  significant 
mixing  in  the  oceanic  boundaries  is  believed  to  be  due  to  internal  wave  reflection  from 
sloping  terrain  [Eriksen,  1985,  1998.  Slinn  and  Riley,  1996.  1998].  The  process  arises  when 
an  oncoming  internal  wave,  traveling  at  some  angle  ft  to  the  horizontal,  reflects  from  a 
sloping  terrain  a  that  is  nearly  equal  to  ft.  In  this  case,  linear  theory  predicts  the  smaller 
amplitude  wave  is  reflected  with  larger  amplitude,  which  can  breakdown  causing  mixing  to 
occur.  This  process  is  illustrated  in  figure  1.1. 

Eriksen  [1982]  has  made  observations  of  internal  wave  fields  near  sloping  topography, 
noting  the  likely  role  of  internal  wave  induced  mixing  in  boundary  layers  of  approximately 
100m.  The  amplitudes  of  reflected  waves  have  been  observed  to  be  much  less  than  linear 
theory  prediction,  which  is  attributed  to  frictional  dissipation  near  the  boundary.  Labo¬ 
ratory  experiments  by  Cacchione  and  Southard  [1974]  and  Cacchione  and  Wunsch  [1974] 
demonstrate  that  amplification  of  reflecting  internal  waves  are  in  agreement  with  linear  the¬ 
ory  for  cases  away  from  the  critical  angle.  When  the  angle  of  reflection  is  near  the  critical 
angle,  the  amplitude  is  much  less  than  predicted.  Ivey  and  Nokes  [1989]  and  Taylor  [1993] 
have  found  that  for  critical  angle  reflection  the  boundary  layer  turbulence  is  either  steady 
or  unsteady  depending  on  the  slope. 

1.2.5  Double  Diffusion 

The  density  of  the  ocean  is  generally  determined  by  two  scalar  parameters,  temperature 
and  salinity,  with  molecular  diffusivities  two  orders  of  magnitude  different  (e.g.,  typical 
temperature  and  salinity  Schmidt  numbers  for  seawater  are  taken  to  be  Sct  =  ( v/Dy ) 
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Figure  1.1:  Diagram  of  internal  gravity  waves  reflecting  from  sloping  terrain.  As  a  reaches 
a  critical  angle,  the  wave  reflection  has  a  larger  amplitude  than  the  incident  wave,  which 
leads  to  instability  and  wave  breakdown.  (Taken  from  Slinn  and  Riley  [1996],  Fig.  1) 


=  7  and  Scs  =  (v/Dt)  =  700.  This  indicates  that  salt  diffuses  much  more  slowly  than 
temperature).  Also,  the  contribution  of  temperature  to  density  is  the  opposite  of  salt. 
That  is,  when  a  fluid  increases  in  temperature,  its  density  decreases ,  whereas  when  salt 
concentration  increases  the  fluid  density  increases.  This  combination  of  large  diffusivity 
difference  and  opposite  contributions  to  density  can  lead  to  a  double-diffusion  instability 
[Gargett,  2003].  The  classic  example  is  an  area  of  high  temperature,  high  salinity  fluid  over 
an  area  of  low  temperature,  low  salinity  fluid  [Ruddick  and  Gargett,  2003,  Schmitt,  1994]. 
This  scenario  is  common  in  the  subtropical  oceans  where  warm  surface  water  evaporates, 
increasing  the  salinity  level  at  the  surface.  Since  the  region  is  top-heavy  in  salt,  a  fluid  parcel 
from  the  top  will  flow  downward.  As  the  parcel  convects  downward  it  will  exchange  heat, 
but  negligible  salt,  with  the  surroundings  due  to  the  large  difference  in  diffusivities.  Thus, 
the  parcel  remains  denser  than  its  surroundings  and  will  continue  to  accelerate  downward. 
Conversely,  a  fluid  parcel  gaining  heat  will  become  less  dense  and  convect  upward.  This 
process  leads  to  the  formation  of  thin  “fingers”  of  salt  transport  seen  in  shadowgraphs 
[Stern,  I960].  Such  salt  fingers  can  be  found  in  the  ocean  at  the  interface  of  layers  in  a 
thermohaline  staircase. 

Thermohaline  Staircase 

A  thermohaline  staircase  is  a  series  of  uniform  temperature  and  salinity  of  layers,  separated 
by  thin  layers  or  “sheets”  of  high  temperature  and  salinity  gradients.  They  are  associated 
with  double-diffusive  instabilities  since  similar  staircases  have  been  found  in  laboratory 
experiments  [Kelley,  1987,  Stern,  1969,  Turner,  1973].  Thermohaline  staircases  have  been 
documented  in  nature  at  several  locations,  including  the  Tyrrhenian  [Molcard  and  Tait, 
1977]  and  Mediterranean  Seas  [Schmitt,  2003],  the  outflow  of  the  Mediterranean  sea  into 
the  Atlantic  [Elliot  and  Tait,  1977,  Williams,  1974]  and  in  the  subtropical  Atlantic  Ocean 
[Boyd,  1989,  Lambert  and  Stnrges,  1977,  Schmitt  et  al.,  1987].  Data  from  the  subtropical 
Atlantic  Ocean  from  1960’s  to  1990’s  suggest  that  the  staircases  are  permanent  features  of 
the  waters  where  they  are  located  [Schmitt,  2003]. 

The  Caribbean-Sheets  and  Layers  Transects  (C-SALT)  program  was  undertaken  in  1985 
to  examine  staircase  formation  and  microstructure  in  the  tropical  North  Atlantic  off  the 
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coast  of  Barbados.  Over  an  area  of  nearly  1  million  km2,  approximately  10  layers  of  uniform 
temperature  and  salinity  are  generally  observed  between  the  depths  of  200-800m,  each  layer 
ranging  from  5-40m  in  depth.  The  interface  between  layers  is  approximately  l-10m  deep 
in  which  very  large  temperature  and  salinity  gradients  exist  [Schmitt  et  ah,  1987].  The 
temperature  difference  in  the  interfaces  between  layers  ranges  between  0.5  —  0.8°C,  and 
salinity  change  across  the  interface  is  typically  0. 1-0.2  psu  (practical  salinity  units).  It  is 
noted  that  these  differences  are  much  larger  than  those  reported  from  the  Tyrrhenian  (0.1°C 
and  0.03  psu)  and  Mediterranean  (0.2°C  and  0.03  psu)  seas. 

Microstructure  measurements  within  the  C-SALT  layer  interfaces  show  several  “subin¬ 
terface”  high  gradient  zones  1-10  cm  thick  [Gregg  and  Sanford,  1987,  Mannorino,  1987, 
Schmitt  et  ah,  1987].  Salt  fingering  occurs  in  these  subinterfaces,  and  is  believed  to  be  the 
dominant  mechanism  of  mixing  and  overturning  of  the  thicker  layers  due  to  the  high  levels  of 
mixing  vs.  the  small  amount  of  turbulence  measured  in  these  regions  [Gregg,  1989].  The  en¬ 
hanced  mixing  provided  by  salt-fingering  has  implications  ranging  from  weather  prediction 
to  nutrient  replenishment. 

The  strength  of  a  double-diffusive  interface  is  given  by  the  density  ratio  Rp: 


(1.9) 


where  Tz.  Sz  are  the  vertical  temperature  and  salt  gradients,  a  is  the  thermal  expansion 
coefficient  defined  in  (1.3),  and  0  is  the  haline  contraction  coefficient: 


5-1 

POS 


r.v 


Favorable  conditions  for  salt  fingering  occur  when 


1  <  Rp  <  «  100 

ns 

where  ky.  ks  are  the  molecular  heat  and  salt  diffusivities  It  is  noted  that  while  this  range 
seems  large,  the  growth  rate  of  salt  fingers  does  not  become  significant  until  Rp  <  2.0,  with 
theoretical  maximum  growth  rate  at  Rp  ^  l.G.  Thus,  while  most  of  the  ocean  is  favorable 
to  salt  fingering  in  the  sense  that  Rp  >  1.0,  regions  with  1.0  <  Rp  <  l.G  are  more  likely  to 
show  staircase  profiles  [Schmitt,  1981,  1988]. 


Differential  Diffusion 

Double  diffusion  is  a  potential  energy  effect.  Differential  diffusion  is  a  kinetic  energy  effect. 
The  underlying  cause  is  same,  that  salt  diffuses  about  100  times  slower  than  temperature. 
However,  in  the  presence  of  turbulence,  the  two  are  often  assumed  to  mix  at  the  same  rate 
(hence  have  the  same  eddy  diffusivities).  This  assumption  is  based  in  the  limit  of  infinite 
Reynolds  number.  However,  much  of  the  ocean  is  mixed  by  patches  of  finite  Reynolds,  finite 
duration  turbulence  [Mourn,  199G,  Smyth  et  ah,  2005].  In  these  areas  of  finite  Reynolds 
number,  it  is  certainly  plausible  that  mixing  of  temperature  and  salt  occurs  at  different 
rates  due  to  the  large  difference  in  diffusivities  between  the  two.  This  difference  in  mixing 
by  turbulence  is  referred  to  as  differential  diffusion  and  has  implication  in  vertical  diffusion 
of  nutrients  across  the  thermocline  as  well  as  the  accuracy  of  current  mixing  models. 
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Turner  [1908]  is  first  credited  with  demonstrating  differential  diffusion  by  mixing  a  fluid 
stratified  by  either  temperature  or  salt  (independently)  via  oscillating  grid  turbulence  and 
measuring  very  different  turbulent  diffusivities.  Altman  and  Gargett  [1990]  found  similar 
results  when  varying  both  temperature  and  salinity  simultaneously.  Jackson  and  Rehmann 
[2003]  found  that  the  ratio  of  salinity  to  thermal  turbulent  diffusivities  Ks/Nt  is  between 
0.5  and  1  for  buoyancy  Reynolds  number  50  <  Reb  <  500,  values  observed  in  the  ocean. 

Numerically,  the  issue  of  differential  diffusion  is  difficult  to  examine  because  of  the  range 
of  scales  required  to  be  resolved  [Gargett  et  ah,  2003].  For  instance,  in  a  stratified  flow  the 
largest  scale  (smallest  wavenumber)  that  needs  to  be  resolved  is  the  Ozmidov  scale  L0.  The 
small  scales  are  dependent  on  the  scalar  diffusivities,  with  the  wavenumber  to  be  resolved 
determined  by  the  Batchelor  wavenumber 


ks  — 


Sc\k 


V 


where  kv  —  (e//)3)1/4  is  the  Kolmogorov  wavenumber.  For  seawater,  to  model  the  diffusion 
of  temperature  would  require  Sep1/2  =  71/2  ^  2.6  times  the  resolution  than  if  the  Schmidt 
number  was  taken  as  1.  The  increase  in  resolution  required  for  salinity  is  26  times  (Scs1//2  = 
700ly/2  ^  26).  This  can  be  a  serious  limitation  in  numerical  simulations,  particularly  direct 
numerical  simulations.  As  a  result,  tradeoffs  must  be  made  between  Reynolds  number  and 
Schmidt  number  depending  on  the  available  computing  resources  (this  is  the  case  for  all 
flows,  not  specific  to  differential  diffusion).  Gargett  et  al.  [2003],  using  numerical  simulations 
with  a  ratio  of  kr/ks  —  0-1  (rather  than  0.01)  and  maximum  Re  =  C/oTo/Sct  =  99,  show 
turbulent  diffusivity  of  temperature  to  be  greater  than  salinity  by  up  to  22%.  Smyth 
et  al.  [2005]  demonstrate  differential  diffusion  in  simulations  of  mixing  in  breaking  Kelvin- 
Helmholtz  billows,  utilizing  a  ratio  kr/ks  =  0.14.  They  note  that  the  turbulent  diffusivity 
is  dependent  on  Reb,  and  the  ratio  of  turbulent  diffusivities  becomes  1  when  Reb  «  0(  102). 


Chapter  2 


Theoretical  Considerations 


2.1  Equations  of  Motion 

Fluid  flow  is  described  by  equations  of  motion  for  mass,  momentum,  and  internal  energy 
and  an  equation  of  state  for  the  fluid.  These  equations  have  the  general  form: 


^  +  V-(pv)  =  0 
at 

(2.1a) 

+  V  •  (pvv)  +  20  X  (pv)  =  V  •  n  +  pA 
i)t 

(2-lb) 

+  V  •  (pve)  =  V  •  (11  •  v)  —  V  •  q  +  pv  •  A 
dt 

(2.1c) 

P  =  p(CuC2,Ci...), 

(2. Id) 

where  p  is  the  density  of  the  fluid,  v  =  (u\ ,  U2-  ^3)  is  the  velocity  vector,  S2  =  (0,  ft  cos(<£), 
$2  si n(0))  is  the  inertial  frame  rotation  rate.  II  is  the  total  stress  tensor,  A  is  an  external 
acceleration  applied  to  the  fluid,  c  —  E+^uf  is  the  sum  of  internal  energy  U  and  mechanical 
energy,  q  is  the  heat  flux,  and  p  is  the  thermodynamic  pressure.  Equation  (2.1b),  with 
the  Coriolis  term  2f2  x  (pv)  present,  is  written  with  a  rotating  frame  of  reference,  as  is 
often  the  case  when  studying  geophysical  flows.  Equation  (2. Id)  is  an  equation  of  state 
related  to  the  nature  of  the  fluid.  Cx  refers  to  quantities  that  determine  the  state  of 
the  density.  For  example,  in  air  Ci,  Ci ,  C3  might  represent  pressure,  temperature,  and 
humidity,  while  in  sea  water  the  Cfs  might  represent  pressure,  temperature,  and  salinity.  In 
their  present  form,  equations  (2.1)  are  a  set  of  coupled  equations  that  completely  describe 
the  fluid  motion.  However,  solutions  to  equations  (2.1)  are  nearly  impossible  to  obtain,  and 
certain  assumptions  must  be  made  to  simplify  the  equations.  Such  assumptions  include 
assuming  a  Newtonian  fluid,  incompressibility,  and  the  Boussinesq  assumption,  each  of 
which  is  explained  below. 

2.1.1  Newtonian  Fluid 

In  order  to  use  equations  (2.1),  something  must  be  known  or  assumed  about  the  molecular 
forces  that  describe  the  total  stress  tensor  IF  Under  most  circumstances,  especially  in 
geophysical  flows,  the  fluid  is  assumed  to  obey  Newton’s  law  of  viscosity,  originally  described 
by  Newton  in  1687.  Newton’s  law  of  viscosity  includes  two  assumptions: 
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1.  II  is  a  linear  function  of  velocity  gradients  and  thermodynamic  state  II(p,  e,  Vv). 

2.  II  is  symmetric. 

The  first  assumption  is  made  because  it  is  desirable  to  have  II  invariant  under  Galilean 
transformation  (i.e.,  II  is  the  same  in  different  inertial  coordinate  systems).  Velocity  gra¬ 
dients  are  invariant  under  Galilean  transformation,  while  the  velocity  vector  itself  v  is  not 
Galilean  invariant.  The  second  assumption  comes  about  by  assuming  that  the  fluid  surface 
cannot  support  a  moment,  which  forces  II  to  be  symmetric. 

To  begin  describing  II  mathematically,  utilize  the  fact  that  pressure  acts  normal  to  the 
surface  of  a  fluid  element.  Thus,  pressure  can  only  reside  on  the  diagonal  of  II,  allowing  IT 
to  be  written  in  the  form 

n  =  -pi  +  r  (2.2) 

where  I  is  the  identity  tensor,  and  f  is  the  viscous  stress  tensor.  Pressure  is  negative 
because  in  this  case  because  it  acts  as  a  compressive  stress.  Since  II  has  been  assumed  to 
be  symmetric,  and  pressure  is  a  normal  force  acting  normal  to  the  fluid  surface  acting  only 
in  the  diagonal  terms,  it  follows  that  f  must  be  symmetric.  Also,  from  early  observations  it 
is  assumed  that  shear  stresses  in  the  fluid  are  linearly  related  to  velocity  gradients.  Utilizing 
this  information,  a  constitutive  equation  for  f  can  be  formed  [Panton,  1996,  p.130]: 

f  =  2fiE  +  A(V  •  v)I  (2.3) 

where  ft  and  A  are  the  first  and  second  coefficients  of  viscosity  respectively,  and  E  is  the 
symmetric  velocity  gradient  tensor 

E=I(Vv  +  (Vv)t),  (2.4) 

with  superscript  T  indicating  transpose.  Equation  (2.3)  is  referred  to  as  the  Newtonian  con¬ 
stitutive  model.  Substituting  (2.3)  into  (2.2)  yields  the  total  stress  tensor  for  a  Newtonian 
fluid 

n  =  -pi  +  2/1E  +  A(V  •  v)I  (2.5) 

Often  (2.5)  is  simplified  simplified  by  using  Stokes’  hypothesis,  which  states  that  thermo¬ 
dynamic  and  mechanical  pressure  are  equal.  This  leads  to  Stokes’s  assumption  of  A  =  —  |/i. 

Now  consider  the  body  force  pA  and  heat  flux  q  within  the  flow.  For  density  stratified 
flows,  typically  the  only  body  force  assumed  to  act  on  the  flow  it  that  due  to  gravity  (i.e., 
centripetal  acceleration  S2xSixfis  negleged),  yielding: 

pA  =  pg,  (2.6) 

while  heat  flux  is  assumed  to  obey  Fourier’s  heat  conduction  law: 

q  =-kS7T  (2.7) 

where  g  is  gravitational  acceleration  (pointing  downward),  k  is  thermal  conductivity,  and  T 
is  temperature.  Substitution  of  (2.5),  (2.6)  (2.7)  into  (2.1)  yields  the  Newtonian  equations 
of  motion  for  stratified  flows 

|?+V-(pv)  =  0  (2.8a) 


— ^  +  V  •  (pvv)  +  212  x  (pv)  =  -Vp  +  2V  •  (pE)  +  V(AV  •  v)  +  pg 

at 


(2.8b) 
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{  +  V  •  {pvt)  =  -(V  •  pv)  +  V  •  (f  •  v)  +  V  •  (fcVT)  +  p(v  •  g)  (2.8c) 

at 

p  =  p(Ca,C2,C3...)  (2.8d) 

where  the  internal  energy  equation  is  left  in  terms  of  the  viscous  stress  tensor.  Equa¬ 
tions  (2.8),  while  telling  us  something  about  the  kinematic  relationship  between  shear  and 
velocity  gradients,  are  still  as  set  of  coupled  equations  for  p,  v,  U.  A  limited  number  of 
exact  solutions  exist,  and  further  simplification  is  required  to  obtain  more  general  solutions. 
The  next  section  addresses  transforming  (2.8c)  to  an  internal  energy  equation  in  terms  of 
temperature. 


2.1.2  Internal  Energy 

The  total  energy  equation  (2.8c)  is  separated  into  kinetic  and  internal  energy  parts.  The 
kinetic  energy  equation  is  obtained  by  taking  the  dot  product  of  (2.8b)  with  v, 


1  D(pv2) 

2  Dt 


— v  •  Vp  +  v  •  ( V  •  t )  +  p(v  •  g) 


(2.9) 


where  TI  —  2L  +  v  •  V  is  the  material  derivative,  and  the  viscous  stress  tensor  has  not  been 
Dt  ot 

expanded  (it  will  be  shown  later  on  that  the  viscous  term  can  be  neglected).  Note  the  loss 
of  the  Coriolis  term.  This  is  because  the  quantity  2 Ft  x  (pv)  is  JL  v,  hence  v  •  [212  x  (pv)] 
is  zero.  Subtracting  (2.9)  from  (2.8c)  gives  the  internal  energy  equation 


DjpU) 

Dt 


—  — pV  •  v  +  f  :  Vv  +  V  •  (kVT) 


(2.10) 


where  the  viscous  terms  have  been  left  as  f  :  Vv  for  now. 

Since  temperature  c  an  be  measured  directly,  it  is  desirable  to  have  an  equation  for 
temperature  T  instead  of  U.  To  do  this,  first  redefine  (2.10)  in  terms  of  enthalpy  II  = 
U  +  (p/p)  [Bird  et  al.,  2002,  p.  337]: 

pSH  =  f:  Vv  +  V  •  (itVT)  +  5?  .  (2.11) 

Dt  Dt 

Assuming  II  is  a  function  of  p  and  T  only  (as  is  the  case  for  a  Newtonian  fluid),  thermo¬ 
dynamic  equilibrium  of  II  reveals  [Panton,  199G,  p.  27] 


dll 


or)  p 


CvdT  + 


P2 
1 


\  dp  J 

OT 


dp 


pj 


dp 


CpdT  +  ~:[\  +  T0}dp 


(2.12) 


with  the  thermal  expansion  coefficient  (3  defined  as 


P 


(2.13) 
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Multiplying  (2.12)  by  p ,  and  equating  the  right  hand  side  of  (2.11)  to  (2.12)  yields  the 
differential  equation  for  temperature  change: 


pCp^T  =  f  :  Vv  +  V  •  (kVT )  +  fp^z  (2.14) 

P  Dt  Dt 

It  can  be  shown  through  scaling  analysis  that  the  viscous  term  is  negligible.  Defining 
velocity  scale  U  and  length  scale  L,  the  ratio  of  the  major  viscous  term  to  inertial  terms: 

=  (215) 

pCp(DT /Dt)  pCpUT/L  C„  TL 


(2.15)  is  typically  C?(10  7)  [Kundu  and  Cohen,  2002,  p.  120],  allowing  the  viscous  term  to 
be  neglected.  The  equations  of  motion  (2.8)  can  then  be  written  as 


dp 

dt 


+  v  '  (pv)  =  0 


(2.16a) 


dpv 

~dt 


+  V  •  (pvv)  +  21  i  x  (pv)  =  — Vp  +  2V  •  (pE)  +  V(AV  •  v)  +  pg 


pCp°J,-  =  V  •  (fcVT )  +  TpSP- 
Dt  Dt 


(2.1Cb) 


(2.16c) 


p  =  p(p,f,5).  (2.16d) 

Treating  the  atmosphere  as  an  ideal  gas  with  equation  of  state  p  =  pRT,  (3  can  be  shown 
to  be  the  reciprocal  of  temperature 


0  = 


1  dp 

Pdf 


Thus,  fJT  =  1,  reducing  (2.16c)  to: 


pcp 


DT 

Dt 


V  •  (A  VT)  +  Sf  . 

Dt 


(2.17) 


(2.18) 


Alternatively,  using  the  relation  Cp  —  Cv  =  It ,  the  equation  of  state  for  an  ideal  gas,  and 
continuity  in  the  form  =  — p(V  •  v),  (2.18)  can  be  written  as  [Bird  et  ah,  2002,  pg.  337]: 


pCvj^  =  V-(kff)-p(V-v). 


(2.19) 


2.1.3  Static  Stability 

A  column  of  fluid  is  said  to  be  stable  when  higher  density  parcels  are  located  below  lower 
density  parcels  (i.e.,  density  decreases  with  height).  Since  density  is  a  function  of  both 
temperature  and  pressure,  compressibility  effects  may  become  important  in  determining 
static  stability.  For  example,  consider  a  fluid  element  at  some  height  z\  that  is  displaced 
adiabatically  (no  heat  transfer  with  the  environment)  to  some  lower  elevation  (and  higher 
pressure)  .  As  the  fluid  element  is  moved  to  the  higher  pressure  region  it  will  become 
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compressed.  The  pressure  has  done  work  on  the  element,  and  the  element’s  temperature 
will  increase.  This  is  shown  mathematically  using  the  first  law  of  thermodynamics: 

dE  =  dQ  +  dw  (2.20) 

where  dE  is  change  in  energy,  dQ  is  change  in  internal  heat,  and  dw  is  work. 

In  order  to  determine  static  stability  of  a  column  of  fluid  a  conserved  quantity  must  be 
used.  A  conserved  quantity  remains  the  same  when  it  undergoes  an  adiabatic  process,  i.e., 
a  process  by  which  dQ  =  0.  To  see  if  temperature  is  conserved,  begin  by  noting  that  the 
work  clone  is  related  to  the  change  in  volume: 

dw  =  —p  ♦  do ,  (2.21) 

where  dv  is  the  change  in  specific  volume  (u  =  4).  Conservation  of  temperature  is  examined 
for  the  atmosphere,  where  is  air  is  assumed  to  obey  the  ideal  gas  equation  of  state  (and 
thus  neglecting  humidity).  The  same  analysis  can  be  performed  for  seawater,  but  with  a 
much  more  complicated  equation  of  state  (see  Appendix).  For  an  ideal  gas,  internal  energy 
is  defined  as 

E  =  Cvf .  (2.22) 

Substituting  (2.21)  and  (2.22)  into  (2.20)  (with  dQ  =  0)  yields: 


CvdT  = 

0  —  p  •  dv 

CvdT  = 

~f>d(  — 
\  P 

CvdT  = 

Rdf  + 

CpdT  = 

RTdp 

P 

dT 

_  dp 

K~ 

RTdp 

P 


T  P 


(2.23) 


where  k  =  Tt/Cp ,  and  the  relationship  Cp  —  Cv  +  R  lias  been  utilized.  Integrating  (2.23) 
from  a  reference  temperature  and  pressure  (usually  taken  at  sea  level)  T0,  pQ  gives: 


f_ 

T0 


(2.24) 


demonstrating  that  temperature  is  not  conserved  when  an  adiabatic  process  is  applied. 

Since  temperature  is  not  a  conserved  quantity,  another  quantity  must  be  used  to  deter¬ 
mine  the  stability  of  the  atmosphere.  Potential  temperature  0  is  conserved,  and  is  defined  as 
the  temperature  a  fluid  particle  would  be  if  it  were  adiabatically  compressed  (or  expanded) 
from  its  in  situ  pressure  p  to  the  reference  pressure  pQ : 


16 


CHAPTER  2.  THEORETICAL  CONSIDERATIONS 


To  see  if  0  is  conserved  under  an  adiabatic  process,  start  by  rearranging  (2.25)  and  differ¬ 
entiating: 


T 

df 


K 


K 


+  K0 


Substituting  these  into  (2.23)  shows  that  potential  temperature  is  conserved  under  an  adi¬ 
abatic  process: 


de 

© 


=  0. 


(2.26) 


Thus,  the  gradient  of  potential  energy  dS/dz  is  used  to  determine  the  stability  of  the 
atmosphere,  with  positive  dQ/dz  signifying  a  stable  stratification. 

Compressibility  effects  are  important  in  the  ocean,  particularly  the  deep  ocean  where 
large  hydrostatic  pressures  exist.  Since 


2.1.4  Incompressible  Flow 


The  assumption  of  incompressibility  does  not  necessarily  mean  constant  density  flow.  Rather, 
the  incompressible  flow  assumption  is  that  changes  in  density  of  a  fluid  particle  are  negligi¬ 
ble.  To  see  when  incompressibility  can  be  applied,  first  assume  that  density  is  a  function  of 
pressure  and  temperature  (neglecting  salinity  and  humidity),  which  yields  a  general  ther¬ 
modynamic  transport  equation  for  density  [Panton,  1996,  p.  230]. 


1  Dp  _  Dp 
P  Dt  Dt  "  Dt  ’ 

where  5  is  the  isothermal  compressibility  coefficient  defined  as 


T\-  1  dP 

a(p,T)  =  -  — 

pop 


(2.27) 


(2.28) 


and  j3  is  the  bulk  or  thermal  expansion  coefficient  defined  in  (2.13). 

The  next  step  is  to  nondimensionalize  equations  (2.27)  and  (2.16).  Denoting  L,  U,  p0 
as  length,  velocity,  and  density  scales,  the  following  nondimensional  quantities  are  defined: 


v 

~  U 


p,U2/k 


_  t, 

~  L/U 


P  = 


P_ 

Po 


Substituting  the  above  scales  into  (2.27)  yields  the  following  nondimensional  density  trans¬ 
port  equation: 


1  Dp  .  r2 

~~FR  =  1M 

p  Dt 


a 


Dp 

~Dt 


PrB  (3DT 
~A  Dt 


(2.29) 


while  substitution  of  the  above  scales  into  (2.16)  yields  the  following  nondimensional  equa¬ 
tions  (where  equation  of  state  is  omitted): 


dp 

dt 


+  V  •  (pv)  =  0 


(2.30a) 
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Opv 


m  +v-(/^v)  +  -(nxv)  =  -vp+- 


2V  •  GuE)  +  V(AV  •  v)  ]  +  pp  (2.30b) 


DT 

~Dt 


pCp  RePr 


k\72T  +  (3B 


1  7  A/ 2  _ 

- h  - - T 

Pr  A 


Dp 

~Dt' 


(2.30c) 


In  the  above  equations,  the  the  following  non-dimensional  quantities  have  been  obtained: 


Re  = 

F2  = 


pLU 

A 


a 


_Cv 

7  Cv 


A  =  apCpT  B  =  6T . 


Here  /  =  2S2sin(0)  is  the  Coriolis  frequency,  and  a  is  the  speed  of  sound  in  the  fluid  medium, 
defined  as 


a  = 


Op 

dp 


s 


From  (2.29),  it  is  evident  that  for  change  in  density  to  be  small,  the  square  of  the  Maeli 
number  M2  must  be  small.  Also,  if  one  considers  a  flow  decelerating  from  some  velocity  v 
to  0,  the  pressure  change  will  be  A p  =  \pv.  Substitution  into  M2  yields  [Panton.  1996,  p. 
236] 


M2  = 


v2 


~2dP 
2  dp 


pv1  Op 
p  Op 


A-1  Ap 

pAp 

Ap 

P 


(2.31) 


Hence,  M2  can  be  seen  as  the  relative  change  of  density  of  the  fluid.  For  incompressibility 
(negligible  density  change).  M2  must  be  small,  and  terms  containing  M2  can  be  neglected. 
In  addition,  nondimensionalizing  thermodynamic  properties  (e.g.,  p,  Cp)  shows  that  when 
the  Mach  number  is  small  relative  changes  in  these  properties  are  also  small,  allowing  the 
properties  to  be  considered  constant.  When  a  low  Mach  number  is  assumed,  terms  with 
M2  can  be  neglected  in  (2.29),  resulting  in  the  following  density  transport  equation: 


(2.32) 


Substituting  (2.33d)  into  (2.30)  and  neglecting  terms  with  M2  yields  the  following  set  of 
(dimensional)  equations  of  motion: 


V  •  v  =  0 


(2.33a) 


p(— -  +  v  •  Vv)  +  2p(El  x  v)  = 
Ot 


vp- t  p  \ 


pCv—  =  kV2T 
Dt 

P  =  P(CUC2,C3...), 


(2.33b) 


(2.33c) 


(2.33d) 


18 


CHAPTER  2.  THEORETICAL  CONSIDERATIONS 


where  (2.33a)  is  obtained  via  substitution  of  (2.32)  into  the  full  continuity  equation  (2.16a), 
(2.33b)  results  from  substitution  of  the  simplified  continuity  equation  (2.33a),  and  (2.33c) 
is  obtained  via  use  of  (2.19)  and  substitution  of  continuity  equation  (2.33a).  Also,  with 
jx  constant  the  viscous  term  simplifies  to  /iV2v  by  substituting  of  equation  2.4  (in  indicial 
notation): 


2/7  V  •  E 


10  f  diii  Diij  \ 
1 2  Oii  \  dxj  ^  dii ) 
__  /  02Ui  <)2  fij  \ 

^  \Oxidxj  ^  dxiOii) 

/7  [V(V  •  v)  +  V2v 
/7V2v 


(2.34) 


where  V(V  •  v)  is  zero  from  continuity. 

The  set  of  equations  (2.33)  are  termed  the  ineompressible  flow  equations  and  are  tyj>- 
ieally  considered  valid  for  M  <  0.3.  A  benefit  of  assuming  ineompressible  flow  (low  Mach 
number)  is  that  the  energy  equation  (2.33c)  is  decoupled  from  the  mass  and  momentum 
equations.  This  is  important,  as  it  allows  for  solution  of  four  variables  (pressure  and  three 
components  of  velocity)  with  four  equations  (continuity  and  three  momentum  equations) 
without  solving  for  temperature  or  internal  energy. 


2.1.5  Boussinesq  Approximation 

The  Boussinesq  approximation  is  a  widely  applied  approximation  to  the  equations  of  motion 
first  suggested  by  Boussinesq  11903].  In  words,  the  Boussinesq  approximation  involves  two 
assumptions: 

1)  Density  fluctuations  within  fluid  motion  are  the  result  of  thermal  effeets  only  (i.e.,  no 
pressure  effeet,  salinity,  humidity). 

2)  Accelerations  within  the  fluid  are  small  compared  to  acceleration  due  to  gravity.  Thus, 
density  fluctuations  are  unimportant  in  the  flow  unless  multiplied  by  gravity. 

Following  the  presentation  by  Spiegel  and  Veronis  [1960],  begin  by  representing  density  in 
the  following  form 

Pi{x,  y,  2,  t)  =  po  +  f>{z )  +  p(x,  y,  z,  t)  (2.35) 

where  po  is  the  constant  spatial  average  density,  p(z)  is  the  variation  of  density  in  the 
absenee  of  motion,  and  p  is  the  fluctuation  of  density  resulting  from  fluid  motion.  The  scale 
height  of  density  is  defined  as 

d  = 

P  po  dz 

In  the  Boussinesq  approximation  it  is  assumed  that  the  fluid  motion  is  limited  to  a  layer  of 
thickness  d  that  is  much  less  than  the  scale  height  Dp  [Spiegel  and  Veronis,  1960],  i.e., 


(2.36) 
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Integrating  (2.37)  from  minimum  to  maximum  density  within  d  yields 

=  e  1  (2.38) 

Po 


where  A p  is  the  maximum  density  variation  in  the  layer  d.  In  addition  to  (2.38),  the  motion 
induced  density  fluctuations  p  are  restricted  to  not  exceed  the  static  variation,  i.e., 


P_ 

Po 


<  0(e)  . 


(2.39) 


Spiegel  and  Veronis  [1960]  note  that  while  (2.39)  must  be  verified  after  solutions  have  been 
obtained,  there  has  been  no  experimental  evidence  that  p  ever  exceeds  A pv.  Using  these 
criterion,  substituting  (2.35)  into  continuity  equation  yeilds: 


V-v 


D 


Dt 

d 


--{po  +  P  +  p) 


P 


P 


1  +  ^  +  4 


Dt 
D_ 

Dt  ' 


Po  Po 


b>+^)+oic!) 


(2.40) 


To  order  c,  (2.40)  can  be  written  as 

V-v  =  0  (2.41) 

which  is  the  same  result  as  assuming  incompressible  flow. 

Now  consider  the  hydrostatic  equation  of  motion.  Expressing  pressure  in  the  same  form 
as  density  in  (2.35)  and  substituting  into  the  vertical  component  of  momentum  gives  the 
following  hydrostatic  equation 

|U-g  (PO  +  P)  (2.42) 

where  the  fluctuating  density  p  is  zero  since  there  is  no  fluid  motion.  Subtracting  (2.42) 
from  the  total  momentum  equation  (2.16b)  with  the  full  expression  of  density  and  pressure 
gives 

+  V  •  (pvv)  +  2 mi  x  v)  =  -Vp  +  2V  •  (/IE)  +  V(AV  •  v)g p .  (2.43) 

at 

Using  the  simplification  of  (2.41),  and  dividing  by  p,  (2.43)  can  be  reduced  to 


SX  +  2 FI  x  v  =  --lv/3  -  ge  ~  +  PV2v  .  (2.44) 

Dt  P  Ap 


where  p  and  0  =  /i/p  have  been  taken  to  be  a  constant.  Substituting  (2.38)  for  e,  the 
Boussinesq  equations  of  motion  can  be  written  in  their  usual  form  as: 


V  •  v  =  0 

^r-r  +  2 il  X  V  =  —  —  Vp/  —  g4-  +  t'V2v 
Dt  Po  Po 

pCvS^  =  kS72T 
Dt 


(2.45a) 

(2.45b) 

(2.45c) 
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P  =  p(Ci ,  C2,  C3  . . . )  (2.45d) 

Often  it  is  necessary  to  have  an  evolution  equation  for  the  density  of  the  fluid.  This  can 
be  achieved  through  the  internal  energy  equation  2.45c  and  the  assumption  that  density 
fluctuations  are  not  a  result  of  pressure  effects  (i.e.,  the  first  assumption  made  at  the 
beginning  of  this  section),  allowing  the  pressure  dependence  in  the  thermodynamic  equation 
of  state  for  density  (2.27)  to  be  eliminated.  In  addition,  since  density  fluctuations  are 
assumed  to  be  small  (cf.  (2.38)),  a  linear  relationship  between  the  density  and  temperature 
fluctuation  can  be  made.  This  allows  substitution  of  Dp  for  DT  in  (2.45c),  yielding  the 
following  Boussinesq  equations: 


V  •  V  =  0 

(2.46a) 

Dv 

+  217  x  v  =  ——\7p'  —  g b-  +  £VJv 

(2.46b) 

bt 

Po  Po 

SI  =  svv 

Dt 

(2.46c) 

where  k  =  k/(pCv )  is  the  mass  diffusivity.  Equations  (2.46)  are  those  commonly  solved 
in  numerical  experiments.  As  with  the  low  Mach  number  approximation,  the  Boussinesq 
approximation  allows  mass  and  momentum  to  be  decoupled  from  the  energy  (now  density) 
equation,  providing  five  variables  (pressure,  three  velocities,  density)  with  five  equations 
(continuity,  three  momentum,  density).  However,  the  assumptions  leading  to  the  Boussi¬ 
nesq  simplifications  should  not  be  overlooked,  namely  that  internal  accelerations  are  small 
compared  to  gravity  and  the  fluid  motion  is  in  a  layer  much  less  than  the  scale  height.  If, 
for  example,  simulations  were  to  be  performed  of  the  entire  thermocline,  fluid  motions  may 
exists  in  a  layer  d  which  is  not  negligible  compared  to  the  thermocline  scale  height,  in  which 
case  the  Boussinesq  approximation  will  not  be  valid. 


2.1.6  Nondimensional  Boussinesq  equations 

Nondimerisionalization  of  the  equations  of  motion  (2.46)  can  provide  insight  to  the  relative 
importance  of  each  term.  In  order  to  nondiminensionalize  (2.46),  the  following  nondimen- 
sional  terms  will  be  used: 


v 

v  =  — 


U 


P 


L 

A l 

Az 

L/U 

P 

PoU2' 


(2.47) 


In  2.47  above,  L  is  a  length  scale,  U  is  a  velocity  length  scale,  and  /  =  212  sin (9)  is  the 
Coriolis  parameter  (where  0  is  the  Earth’s  latitude).  Note  that  in  density  stratified  flows 
it  is  common  to  nondimensionalize  p  by  by  the  product  of  the  length  scale  and  the  scale  of 
the  change  in  density  Ap/Az.  Also  note  that  pressure  is  nondimensinalized  by  the  dynamic 
pressure  poU2.  Substituting  (2.47)  into  equations  (2.46)  yields  the  following  nondimensional 
Boussinesq  equations  of  motion: 

V  •  v  =  0  (2.48a) 

15i  +  l£(n  x  v)  =  -Vp  “  (f )  '*■ +  I^v"v  (2'48b) 
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=  JL  vV 

Dt  ScRe  ' 


(2.48c) 


The  following  nonclimensional  parameters  are  used  to  determine  the  significance  of  each 
terms  in  equatoins  (2.48): 


v 


K 


(2.49) 


where  Ro  is  the  Rossby  number,  F  is  the  Froude  number,  Re  is  the  Reynolds  number,  Sc  is 
the  Schmidt  number,  and  N2  =  (g/ p)Ap/Az  is  the  Brunt-Vaisala  (or  buoyancy)  frequency. 
The  significance  of  each  nondimensional  parameter  is  briefly  described  below. 

The  Rossby  number  (Ro)  describes  the  ratio  of  inertial  to  Coriolis  forces.  When  Ro 
is  large,  Coriolis  forces  due  to  planetary  motion  can  be  neglected.  This  will  occur  in  low- 
latitude  regions  of  the  Earth  (i.e.,  in  the  Tropics),  when  the  length  scale  is  small,  or  when 
the  velocity  is  very  large. 

The  Fronde  number  (F)  describes  the  ratio  of  inertia  forces  to  gravity  forces.  In  a 
stratified  fluid,  a  buoyancy  force  arises  due  to  the  motion  of  fluid  elements  with  differing 
density.  For  instance,  a  fluid  element  with  a  high  density  traveling  upward  to  a  lower  density 
region  would  experience  a  buoyancy  force  trying  to  push  it  back  to  its  lower,  stable  location. 
Decreasing  F  signifies  increasing  stratification  strength.  Note  that  buoyancy  force  can  act 
as  a  stabilizing  force  that  inhibits  vertical  motion. 

The  Reynolds  number  (Re)  describes  the  ratio  of  inertia  forces  to  viscous  forces.  When 
Ro  is  large,  the  inertial  terms  of  the  equations  of  motion  dominate  and  viscous  terms  can 
be  neglected,  such  as  flows  far  from  boundaries  or  with  very  large  length  scale.  Conversely, 
a  small  Re  indicates  a  viscous  dominated  flow. 

The  Schmidt  number  (Sc)  is  used  to  describe  the  ratio  of  momentum  diffusivity  (i.e.. 
viscosity)  to  mass  diffusivity.  A  large  Sc  indicates  that  momentum  is  diffused  faster  than 
mass.  Note  the  Sc  ~  0.7  for  air,  while  in  sea  water  Sc  is  approximately  7  for  temperature 
and  700  for  salinity.  Often  numerical  simulations  are  performed  with  Sc  =  1,  since  a  higher 
resolution  simultation  is  required  to  capture  effects  of  both  viscosity  and  diffusivity  when 
Sc  is  much  different  than  1. 


2.1.7  Low  Froude  Number  Equations 

In  density  stratified  flow,  the  buoyancy  force  of  the  stratification  acts  to  suppress  (positive) 
vertical  motions.  This  leads  one  to  believe  that  the  horizontal  and  vertical  scales  will  be 
different.  As  such,  it  is  common  in  the  literature  to  find  equations  (2.4G)  and  (2.48)  divided 
into  horizontal  and  vertical  components  [Billant  and  Chomaz,  2000a,  Gargett,  1988,  Lilly, 
1983,  Riley  and  Lelong,  2000,  Riley  et  al.,  1981,  e.g].  They  are  usually  written  as: 


~  -  du; 

V//-VH+  ^7  =  0 
az 


Oyfj_ 

Of. 


+  V//  •  Vv//  -j-  w—  +  ft 

az 


Ow  _  ~  _  Ow  Op  _  ~  o  _ 

+  VH  *  vk;  +  =  -Tp  -  P9  +  V  w 

dt.  az  az 


°P  ,  ~  fy  -  -  UP  - 

— ^  +  v//  •  Vp  +  w—  -w  =  V 

Of  dz 


P 


(2.50a) 


(2.50b) 


(2.50c) 


(2.50d) 
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whore  subscript  H  indicates  horizontal  (e.g.,  v//  =  (6,5)),  w  is  vertical  velocity,  all  other 
term  are  as  defined  above. 

Nondimensionalization  of  the  equations  of  motion  provides  insight  to  the  relative  impor¬ 
tance  of  each  term.  Nondimensional  versions  of  equations  (2.50)  are  presented  by  several 
authors,  including  Billant  and  Chomaz  [2000a],  Lilly  [1983],  Riley  and  belong  [2000],  Riley 
et  al.  [1981].  As  mentioned  above,  stratification  typically  suppresses  vertical  motion,  which 
leads  to  different  vertical  and  horizontal  scales.  Following  the  scaling  arguments  of  Riley 
et  al.  [1981],  choose  the  following  length,  horizontal  velocity,  and  time  scales: 


xu  ~  Lh 

Z  ~  Ly 

v//  ~  U 

a  — 

-  lU 72 

W  —  F 

f>  ~  PoU 

Lh 

a 

t  ~  Ljj/u 

_  ,  U 2 

P  ~  PoZT~ 

pLy 


with  Froude  number  F  =  U/(NLy).  The  resulting  nondimensional  continuity,  horizontal 
momentum,  vertical  momentum,  and  density  equations  of  motion  then  become  [Riley  and 
Lelong,  2000]: 

V//  •  uH  +  F2^  =  0  (2.51a) 

az 

tt.vh  +  vw  •  V//v//  +  F2uAuh  +  Ae2  x  v//  =  IIP  +  *4tt-v2v//  (2.51b) 

Ot  az  Ro  o2  Re 


o  o  /  dw  „  9  dw 

a2 F2  (  — — h  v//  •  Vu;  +  F2w-— 

dt  (JZ 


0P  F'2  n2 


(2.51c) 


Op 

Ot 


+  v//  •  Vp  +  F 2w^~  -  w  = 
az 


1  1 

o2  ScRe 


V2p 


(2.51d) 


where  Re  =  UL/u  is  a  Reynolds  number  based  on  energy  containing  length  L,  Sc  =  v/k 
is  the  Schmidt  number,  and  Ro  =  UfL  is  the  Rossby  number.  As  stratification  becomes 
important  (F  <  0(1)),  all  the  terms  with  vertical  velocity  w  diminish  in  importance  and  can 
be  neglected.  Evidence  of  this  diminished  vertical  scale  of  motion  is  seen  in  the  formation 
of  horizontal  “pancake  eddies”  in  numerous  numerical  and  experimental  studies  (examples 
are  given  in  the  sections  on  turbulence  in  wakes  and  grid  turbulence).  Equations  (2.51)  are 
often  referred  to  as  the  low  Froude  or  stratified  turbulence  equations. 


2.1.8  Length  Scales  in  Stratified  Flows 


The  smallest  scale  of  turbulent  motion  is  taken  as  the  familiar  Kolmogorov  length  scale 


Lk  = 


(2.52) 


where  v  is  kinematic  viscosity  and  ?  is  kinetic  energy  dissipation  rate.  is  taken  as  the 
smallest  size  of  turbulent  motion  before  being  dissipated  into  internal  energy  (heat).  At  this 
small  scale  the  length  is  typically  assumed  to  be  isotropic  regardless  of  density  stratification, 
although  isotropy  is  disputed  by  Smyth  and  Mourn  [2000b]. 

Buoyancy  tends  to  affect’ the  larger  vertical  scales  of  motion.  The  largest  scale  of  turbu¬ 
lent  motion  before  buoyancy  inhibits  vertical  motion  (or  the  largest  scale  at  which  turbulent 
overturning  can  occur)  is  the  Ozmidov  or  overturning  scale: 


L0 


(2.53) 
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where  N  is  the  Brunt-Vaisala  (buoyancy)  frequency: 


N  = 


9  Op 
Po  Oz  ' 


L0  is  derived  in  [Ozmidov,  19G5]  by  finding  the  length  such  that  buoyancy  and  inertial  scales 
are  equal.  L0  is  typically  applied  to  flows  with  no  mean  shear. 

If  a  flow  is  subject  to  a  mean  shear,  the  Corrsin  scale  is  the  largest  scale  at  which  the 
flow  is  unaffected  by  the  shear,  defined  as: 


Lc 


(2.54) 


where  S  is  a  measure  of  the  applied  shear. 

Another  measure  of  the  overturning  length  is  the  Thorpe  scale,  Lt  (Thorpe,  1977).  It  is 
computed  by  reordering  the  fluid  elements  such  that  the  vertical  density  profile  is  statically 
stable.  is  then  taken  to  be  the  rms  of  the  distance  d  each  element  was  moved  to  make 
the  density  profile  stable. 

L  t  =  (d2)*=drms  (2.55) 

A  measure  of  the  distance  a  fluid  particle  could  be  displaced  if  all  its  vertical  kinetic 
energy  were  converted  to  potential  energy  is  the  buoyancy  scale  L|>: 

Lb  =  (2. 50) 

where  wrm3  is  the  mis  vertical  velocity.  Lb  is  shown  to  be  closely  proportional  by  Smyth 
and  Mourn  [2000a],  and  is  said  to  be  an  upper  limit  of  Lt  [Mourn,  1996]. 


2.1.9  Additional  Nondimensional  Parameters  in  Stratified  Flows 


The  Richardson  number  is  commonly  defined  in  stratified  flows.  Two  common  Richardson 
number  definitions  exist.  The  bulk  Richardson  number  is  simply  the  inverse  F2,  and  is  a 
description  of  the  overall  influence  of  gravity  on  the  flow.  Tin1  local  Richardson  number  is 
more  commonly  used,  and  is  defined  as: 


Ril 


N2 

(dU/dz) 


Ril  is  thus  a  ratio  of  buoyancy  force  to  shear  forces  in  the  flow.  As  Rij  decreases,  the  flow 
is  subject  to  shear  instabilities.  A  common  criterion  used  is  Rij  <  *  for  Kelvin-Helmholtz 
instability. 

The  buoyancy  Reynolds  number  is  defined  by  taking  the  ratio  of  L()  to  Lk  [Gibson,  1980, 
Gregg,  1987,  Smyth  and  Mourn,  2000a]: 


(2.57) 


L0  is  the  maximum  length  at  which  eddies  can  vertically  overturn  before  being  affected  by 
buoyancy  [Ozmidov,  1965,  Smyth  and  Mourn,  2000a],  while  Lk  is  the  maximum  length  scale 
that  occurs  before  energy  is  dissipated  via  viscous  effects.  Reb  can  be  seen  as  a  measure  of 
the  “bandwidth  of  length  scales  available  to  turbulence”  [Gregg,  1987].  When  Reb  ~  0(  1), 
L0  %  Lk  and  vertical  overturning  (3D  turbulence)  does  not  occur  [Riley  and  Lelong,  2000]. 
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Chapter  3 

Taylor-Green  Simulations 


This  chapter  contains  results  of  numerical  simulations  initialized  with  Taylor-Green  vortices 
and  a  linear  density  profile  (i.e.,  uniform  vertical  density  gradient).  An  overview  of  the 
Tavlor-Green  vortices  is  first  given,  followed  by  descriptions  of  the  equations  of  motion  and 
numerical  method.  Results  are  then  presented. 


3.1  Overview 


Direct  numerical  simulations  (DNS’s)  of  flows  initialized  with  Taylor-Green  vortices  were 
analyzed.  The  stratification  of  the  flow  is  constant  in  time,  and  there  is  no  ambient  shear 
(i.e.,  decaying  turbulent  flow).  Taylor-Green  vortices  were  chosen  because  they  can  be 
considered  idealizations  of  flows  resulting  from  laboratory  experiments  in  which  a  rake  is 
pulled  through  a  continuously  stratified  tank  [e.g.,  Fincharn  et  ah.  1996,  Fraud  et  ah.  2005]. 
Note,  however,  that  our  simulations  are  conducted  with  Schmidt  number  =  1.  vs.  Sct  =  7 
for  water.  Briefly,  the  initial  condition  consisted  of  Taylor-Green  vortices  plus  broad-banded 
noise  with  a  level  approximately  10%  of  the  Taylor-Green  vortex  energy.  The  Taylor-Green 
vortices  satisfied  the  following  mathematical  form: 

\rtg  =  U  cos (kz)  [cos(/ci:)  sin (ky),  —  sin(kx)  co s(kjj),  0] 


where  U  determines  the  initial  velocity  scale,  and  L  =  l/h  determines  the  length  scale  for 
this  field.  Thus,  for  all  simulation  results,  velocities  are  nondimensionalized  by  ZV,  lengths 
by  L,  and  time  by  L/U.  The  Froude  number  and  Reynolds  number  characterizing  the 
simulations  are  defined  respectively  as 


Fr,  = 


2nU 
N  L  ’ 


R  UL 
Rei,  =  —  . 


(3.1) 


Simulations  were  run  with  nominal  Froude  numbers  (Fl)  of  2  and  47  and  nominal 
Reynolds  number  (Ren)  between  200  and  9600.  A  summary  of  the  different  cases  is  given 
in  Table  3.1.  As  described  in  Riley  and  de  Bruyn  Kops  [2003],  the  simulated  flow  exhibits 
many  characteristics  of  stratified  turbulence.  The  horizontal  length  scales  grow,  and  the 
vertical  length  scales  decrease  in  time.  This,  combined  with  decoupling  of  the  horizontal 
motions  in  the  vertical  direction,  leads  to  the  formation  of  horizontal  “pancake”  vortices. 
This  behavior  is  evident  after  a  dimensionless  time,  £,  of  about  15.  In  this  study  the  focus 
is  on  time  t  =  20,  when  the  two-dimensional  character  of  the  flow  is  strong  but  the  flow  is 
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Notation 

Fl 

ReL 

Sc 

Nx 

Ny 

Nz 

F2R2 

2 

200 

i 

256 

256 

128 

F2R4 

2 

400 

i 

256 

256 

128 

F2R8 

2 

800 

i 

256 

256 

128 

F2R16 

2 

1000 

i 

256 

256 

256 

F2R32 

2 

3200 

i 

512 

512 

256 

F2R64 

2 

6400 

i 

768 

768 

384 

F2R96 

2 

9600 

i 

1024 

1024 

512 

F4R2 

4 

200 

i 

256 

256 

128 

F4R4 

4 

400 

i 

256 

256 

128 

F4R8 

4 

800 

i 

256 

256 

256 

F4R16 

4 

1600 

i 

256 

256 

256 

F4R32 

4 

3200 

i 

512 

512 

256 

F4R64 

4 

6400 

i 

768 

768 

384 

F4R9G 

4 

9600 

i 

1024 

1024 

512 

Table  3.1:  Conditions  for  simulations  of  quasi-horizontal  vortices.  Nx,  Ny,  and  Nz  are  the 
number  of  grid  points  in  each  direction. 


also  still  very  energetic.  In  particular,  the  relationship  between  kinetic  energy  dissipation 
rate  and  vertical  shear,  the  definition  of  horizontal  length  scale,  and  parameterization  of 
turbulence  are  examined. 

The  Taylor-Green  initial  condition  and  the  periodic  boundary  conditions  result  in  either 
two  or  four  planes  (depending  on  the  size  of  the  numerical  domain)  of  maximum  shear  and 
the  same  number  of  planes  of  minimum  shear.  This  is  different  from  laboratory  experiments 
or  in  natural  settings  (e.g.,  atmosphere  and  ocean)  in  which  the  number  and  spacing  of 
planes  of  high  and  low  shear  adjusts  to  the  flow  conditions.  As  a  consequence  of  this 
prescribed  spacing,  volume  averages  of  shear  and  dissipation  rate  are  difficult  to  interpret. 
Here  we  consider  only  the  planes  of  maximum  shear  (either  two  of  four)  and  denote  the 
average  over  these  planes  by  an  overbar. 


3.2  Equations  of  Motion 

The  flow  fields  are  assumed  to  satisfy  the  incompressible  continuity  and  Navier-Stokes 
equations  subject  to  the  Boussinesq  approximation.  Also,  for  better  understanding  of  the 
underlying  physics,  the  simulations  are  performed  in  a  non-rotating  frame  of  reference  (i.e., 
the  Coriolis  term  has  been  neglected  since  the  simulation  length  scale  is  small  compared  to 
the  scales  acted  upon  by  Coriolis  forces.)  Thus,  the  governing  equations  are  those  of  (2.46) 
with  no  Coriolis  term.  In  nondimensional  form  these  equations  are  written  as:  (c.f.  (2.46)): 


V  •  v  =  0 


<9v  _  S/p 

—  +  v  •  Vv  = - i- 

Ot  p0 


2tt 

Fl 


Pe2  + 


R('l 


V2v 


l+vV"+"'l  =  &kv2^ 


(3.2a) 

(3.2b) 


(3.2c) 
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Here  v  =  ( it ,  v ,  w)  is  the  velocity  vector,  p  and  p  are  the  density  and  pressure  deviations 
from  their  ambient  values,  ez  is  a  unit  vector  in  the  vertical  (z)  direction,  S c.  —  D/D  is 
the  Schmidt  number  (where  D  is  the  mass  diffusivity),  and  Fl  and  Ren  are  Reynolds  and 
Froude  numbers  as  defined  in  (3.1). 


3.3  Numerical  Method 

The  equations  of  motion  are  solved  using  a  pseudo-spectral  technique.  The  boundaries  are 
taken  to  be  periodic,  allowing  trigonometric,  evenly  spaced  interpolation  points  and  easy 
calculation  of  spatial  derivatives  via  Fast  Fourier  Transform  (FFT)  [Trefethen,  2000].  The 
equations  were  advanced  in  time  using  a  third-order  Adams-Bashforth  scheme  with  pressure 
projection  (fractional  step)  method.  Briefly,  the  pressure  projection  method  is  a  two  step 
method  to  solve  for  the  velocity  field.  First,  the  velocity  fields  are  advanced  in  time  without 
the  pressure  term.  Then,  with  the  condition  of  incompressibility,  pressure  is  expressed  as  a 
Poisson  equation  in  terms  of  the  newly  calculated  velocity  fields.  The  Poisson  equation  is 
then  solved  for  pressure,  then  the  value  obtained  for  pressure  is  used  to  modify  or  “correct” 
the  velocity  fields.  Further  discussion  of  the  pressure  projection  method  can  be  found  in 
Perot  [1993]  and  Chang  et  al.  [2002]. 

In  order  to  eliminate  the  majority  of  aliasing  errors,  a  spherical  wave-number  truncation 
of  approximately  15/16  with  Kmax  the  maximum  wave  number  in  the  discrete  Fourier 

transforms,  was  used.  The  momentum  equation  was  advanced  in  time  with  the  nonlinear 
term  expressed  in  vorticity  form,  while  the  alternating  time-step  scheme  suggested  by  Kerr 
[1985]  was  employed  for  the  density  field  to  approximate  the  skew-symmetric  form  of  the 
non-linear  term  and  thereby  minimize  aliasing.  A  skew-symmetric  matrix  is  one  in  which 
the  eigenvalues  are  all  pure  imaginary.  This  implies  the  non-linear  term  advects  without 
causing  growth  or  decay,  and  will  not  cause  the  model  to  fail.  [Boyd,  2001,  pg.213] 


3.4  Simulation  Results 

3.4.1  Relation  between  vertical  shear  rate  and  kinetic  energy  dissipation 
rate 

The  vertical  velocity  in  case  F2R32  at  four  different  times  is  shown  in  Fig.  3.1.  The  white 
bar  in  each  figure  connects  two  material  points  that  are  tracked  in  time.  At  t  =  17.5  the 
flow  is  fairly  quiescent,  but  by  t  =  20  a  turbulent  patch  has  begun  to  form  in  the  vicinity 
of  the  white  line.  This  patch  continues  to  develop  through  time  22.5. 

By  slicing  the  domain  vertically  through  the  white  line,  the  dynamics  of  this  particular 
instability  can  be  studied.  In  Fig.  3.2,  the  horizontal  velocity  in  the  direction  of  the  white 
line  is  plotted  for  each  of  the  four  times.  The  figures  are  colored  so  that  black  indicates 
flow  to  the  left  and  white  indicates  flow  to  the  right.  The  gray  bars  above  each  panel 
in  Fig.  3.2  correspond  to  the  white  lines  in  Fig.  3.1.  At  t  —  17.5  there  is  little  shearing 
action  in  the  region  of  the  gray  bar,  but  by  t  =  20  moderate  shearing  has  developed.  By 
t  —  21.5  the  shearing  is  strong,  and  Kelvin-Helmholtz  roll-ups  are  apparent.  The  roll-ups 
are  even  more  apparent  in  the  corresponding  slice  through  the  total  density  field  at  t  =  21.5, 
published  in  de  Bruyn  Kops  et  al.  [2003].  This  type  of  qualitative  analysis  was  repeated  for 
a  number  of  the  simulation  cases.  While  it  was  not  always  possible  to  find  a  clean  Kelvin- 
Helmholtz  roll-up  associated  with  each  turbulent  patch,  high  vertical  shear  accompanied  the 


28 


CHAPTER  3.  TAYLOR-GREEN  SIMULATIONS 


turbulence  in  all  cases  observed.  Although  this  analysis  does  not  provide  a  definitive  answer, 
it  strongly  suggests  vertical  shear  as  the  dominant  mechanism  for  triggering  turbulence  in 
these  simulations. 

When  developing  theories  and  models  for  turbulence  subject  to  strong  stable  stratifica¬ 
tion,  it  is  common  to  assume  that  vertical  shear  of  the  horizontal  motions  causes  most  of 
the  dissipation  rate  of  kinetic  energy,  i.e., 


uS2  «  e  .  (3.3) 

Here  e  =  2veijCi3  is  the  kinetic  energy  dissipation  rate,  is  the  symmetric  part  of  the 
viscous  stress  tensor  r^,  and  S 2  =  (< du/dz )2  +  ( dv/dz )2  is  the  square  of  vertical  shear.  For 
instance,  Billant  and  Chomaz  [2001]  and  Riley  and  de  Bruyn  Kops  [2003]  make  assumption 
(3.3)  in  order  to  estimate  a  vertical  length  scale;  Shih  et  al.  [2005]  make  this  assumption  in 
order  to  relate  the  buoyancy  Reynolds  number  Reb,  to  the  ratio  of  Reynolds  and  Richardson 
numbers.  Sup>port  for  (3.3)  comes  from  numerical  simulations,  [e.g.,  Herring  and  Metais, 
1989],  and  laboratory  experiments  that  span  a  wide  range  of  Reynolds  numbers.  In  par¬ 
ticular,  the  experiments  of  Finchain  et  al.  [1996]  show  that  vS2 je  «  0.9,  a  result  that  is 
verified  experimentally  by  Fraud  et  al.  [2005]. 

Analysis  will  begin  by  considering  v  ( S 2)  /  (e)  for  each  of  the  cases  in  Table  3.1,  shown 
in  Figure  3.3,  and  each  component  (e)7J  as  a  fraction  of  e,  shown  in  Figures  3.4(a)  and  (b). 
Here  (♦)  denotes  a  volume  averaged  quantity.  For  ReL  =  800  at  both  F[  =  2  and  4,  the 
results  are  entirely  consistent  with  those  of  Fincharn  et  al.  [1996]  and  Fraud  et  al.  [2005], 
in  that  vertical  shear  accounts  for  about  90%  of  the  dissipation  rate.  The  horizontal  shear 
terms  are  very  small,  as  are  the  contributions  to  e  of  the  normal  strain  rates.  There  is  some 
dissipation  due  to  vertical  motion,  particularly  in  the  Fl  =  4  case,  but  it  is  small. 

For  Ren  =  200,  the  flow  conditions  are  markedly  different  with  v  (S2)  /  (e)  ~  0.6.  The 
contributions  of  each  component  of  (e)  suggest  a  two-dimensional  Stokes-like  flow  in  which 
vertical  motion  is  almost  completely  inhibited  by  gravity  for  both  Fronde  number  cases  and 
the  isotropnc  character  of  the  pressure  force  causes  all  of  the  horizontal  contributions  of  (e) 
to  approach  their  isotropic  values  (i.e.,  the  flow  never  becomes  turbulent).  Note  that  in 
isotropic  turbulence,  the  contributions  to  (e)  of  the  diagonal  terms  in  the  strain  rate  tensor 
are  13.3%,  while  those  of  the  off-diagonal  terms  are  10%. 

For  Reynolds  numbers  above  800,  the  relative  contributions  of  each  component  of  (e) 
change  rapidly  with  increasing  ReL-  Fhysical  reasoning  suggests  that  all  components  of 
(e)  that  depend  on  vertical  velocity  will  increase  as  the  flow  becomes  more  turbulent  and 
the  vertical  velocity  increases.  Furthermore,  such  reasoning  suggests  that  this  phenomenon 
will  occur  at  lower  ReL  for  the  Fl  =  4  cases  than  for  the  Fl  =  2  cases  since  the  vertical 
motions  are  more  strongly  suppressed  by  gravity  at  lower  Fronde  numbers.  The  simulation 
data  is  consistent  with  this  reasoning  in  that  (e)33  /  (e)  increases  with  Reynolds  number 
for  both  Fl  =  2  and  4  and  is  higher  for  the  Fl  =  4  cases  except  at  the  highest  Reynolds 
numbers.  Considering  this  result  alone,  we  might  conclude  that  the  simulated  flow  is  insuf¬ 
ficiently  stratified  for  (3.3)  to  hold  for  the  higher  Reynolds  number  cases.  Note,  however, 
the  curves  for  the  (<s)n  and  {s)22  *n  Figures  3.4(a)  and  (b).  The  significant  rise  of  these  two 
contributions  to  (e)  shows  that  it  is  not  just  increasing  vertical  motion  that  causes  vS2 /  (e) 
to  decrease  to  about  0.4  at  Rcl  =  9600.  Dissipation  due  to  normal  strains  contributes 
significantly  to  the  total  dissipation  rate  at  higher  Reynolds  numbers. 

While  it  would  be  convenient  if  v  (S2)  /  (e)  &  1  for  all  cases  of  strongly  stratified  turbu¬ 
lence,  for  most  modeling  and  theoretical  applications  introduction  of  an  order  one  constant 
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t=17.5  t=20 


Figure  3.1:  A  horizontal  slice  through  the  vertical  velocity  field  at  the  plane  of  maxi¬ 
mum  shear  at  four  different  times  in  a  Taylor-Green  simulation  with  Froude  number  2  and 
Reynolds  number  3200.  The  white  bar  connects  two  material  points  that  move  with  time. 
The  points  start  in  a  region  of  relative  calm,  experience  an  instability,  and  end  up  in  a 
turbulent  patch. 
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t=1 7.5 


t=20 


t=21 .5 


t=22.5 


Figure  3.2:  The  horizontal  speed  ori  vertical  planes  aligned  with  the  white  bar  in  Fig.  3.1. 
Black  indicates  negative  and  white  indicates  positive.  The  gray  bar  above  each  panel 
corresponds  to  the  white  bar  in  Fig.  3.1. 


Figure  3.3:  Ratio  of  (v  (S2))  /  (e)  vs.  Rcl-  Fl  =  2  (•)  :  Fl  =  4  (■) 
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Figure  3.4:  Contribution  of  the  six  independent  terms  of  (s)^  normalized  by  total  (e)  vs. 
R('i  for  (a)  F\^2  and  (b)  Fl=4.  The  horizontal  dashed  lines  mark  the  theoretical  values 
for  the  normal  and  shear  components  in  isotropic  turbulence. 


would  be  acceptable  provided  that  S2  and  e  were  well  correlated.  This  correlation  is  con¬ 
sidered  in  Figure  3.5  in  which  local  values  on  the  planes  of  high  shear  of  vS 2  are  plotted 
versus  e  for  two  different  simulation  cases.  In  the  top  panel  of  the  figure,  it  is  apparent  that 
for  ReL  —  800,  not  only  are  the  square  of  the  shear  and  the  dissipation  rate  well  correlated 
but  there  are  very  few  points  far  from  the  diagonal.  For  this  case,  relation  (3.3)  is  not  only 
excellent  on  average,  it  is  excellent  locally  in  regions  of  high  shear.  In  the  bottom  panel 
of  the  figure,  the  results  for  Rei,  =  0400  show  that  relation  (3.3)  is  not  very  good  even  to 
within  a  multiplicative  constant  for  this  case. 

3.4.2  Buoyancy  Reynolds  Number 

In  order  to  understand  why  the  laboratory  experiments  conducted  over  a  wide  range  of 
Reynolds  numbers  consistently  show  v(S2)/e  ~  0.9  while  our  simulations  support  this 
relationship  only  for  a  fairly  narrow  range  of  Reynolds  numbers,  we  consider  now  both 
laboratory  and  simulation  data  in  terms  of  the  buoyancy  Reynolds  number  (2.57).  This 
quantity  has  been  used  extensively  in  the  parameterization  of  stratified  turbulence,  [e.g., 
Gibson,  1980,  Gregg,  1987,  Imberger  and  Boashash,  1986,  Smyth  and  Mourn,  2000a]  and 
can  be  derived  from  the  ratio  of  the  Ozmidov  scale  (2.53)  and  the  Kolmogorov  scale  (2.52). 

For  each  of  the  simulation  cases,  the  planar  average  buoyancy  Reynolds  number,  (Reb)  = 
(e)  /vN2,  is  plotted  versus  Ren  in  the  top  panel  of  Figure  3.6.  Since,  as  shown  by  Riley 
and  de  Bruyn  Kops  [2003],  the  dissipation  rates  for  all  the  cases  are  about  the  same,  (Reb) 
is  very  nearly  proportional  to  Rec-  Note  that  (Reb)  is  computed  at  t  —  20,  whereas  Rec 
is  a  nominal  value  for  the  simulation.  In  the  bottom  panel  of  the  same  figure,  v  (S2)  /  (e) 
is  plotted  versus  (Reb)  for  all  the  simulation  cases.  When  (Reb)  >  1,  the  data  for  both  F^ 
cases  collapse  very  well  onto  a  common  curve  that  decreases  rapidly  with  increasing  (Reb). 
In  the  range  0.1  <  (Reb)  <  1,  there  is  more  scatter  in  the  data  but  high  values  of  v  (S2)  /  (e) 
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are  observed,  consistent  with  results  of  other  numerical  simulations  [e.g.,  Smyth  and  Mourn, 
2000b]  a s  well  as  the  laboratory  results  of  Fine  ham  et  al.  [199G]  and  Fraud  et  al.  [2005]. 

For  the  case  of  Fiiichain  et  al.  [1996]  with  Rejvi  =  6100  and  N  =  2.3  rad  s_1  used  in 
their  Figure  8  to  show  the  contribution  of  vertical  shear  to  dissipation  rate,  the  buoyancy 
Reynolds  number  is  estimated  to  be  about  0.2  at  early  time.  For  the  case  of  Fraud  et  al. 
[2005]  with  Re\i  =  9000  and  Fr\i  =  0.09  used  in  their  Figure  25,  the  buoyancy  Reynolds 
number  is  estimated  to  be  about  0.1  at  early  time.  These  values  of  Re^  correspond  roughly 
to  where  v  (52)  /  (e)  is  maximum  in  Figures  3.6  (a)  and  (b). 

In  both  laboratory  experiments,  however,  e  decreases  by  several  orders  of  magnitude 
over  the  duration  of  the  experiment  with  a  corresponding  decrease  in  buoyancy  Reynolds 
number.  Neither  Fincham  et  al.  nor  Fraud  et  al.  report  a  decrease  in  vS2 /e  as  the 
experiments  evolved.  This  suggests  that  the  behavior  at  low  (Reb)  observed  in  Figures  3.6 
(a)  and  (b)  is  due  to  the  simulated  flows  being  laminar.  Based  on  other  statistics,  we  know 
this  to  be  the  case.  It  is  likely,  based  on  physical  reasoning  and  the  two  sets  of  laboratory 
data,  that  vertical  shear  will  account  for  90%  of  the  dissipation  rate  when  (Reb)  is  order 
one  or  less  provided  that  the  flows  are  turbulent  (of  course,  with  (Reb)  <  1,  turbulence  is 
in  the  the  quasi-2d  sense). 

3.4.3  Horizontal  Length  Scale 

Riley  and  de  Bruyn  Kops  [2003]  postulate  a  length  scale,  L^,  for  horizontal  motions  in  order 
to  define  a  horizontal  Reynolds  number,  R eh,  and  a  horizontal  Froude  number,  F/*,  for  their 
F^Re/j  scaling.  How  L ^  is  to  be  computed  for  practical  use  is  not  discussed  in  that  paper, 
but  their  derivation  of  the  F^Re^  scaling  implies  that  L ^  can  be  defined  in  terms  of  and 
e  (i.e.,  Lh  is  an  advective  length  scale).  Similarly,  when  Reb  is  related  to  the  square  of  a 
Froude  number  and  a  Reynolds  number  [e.g.,  Ivey  and  Imberger,  1991,  Shih  et  al.,  2005],  an 
advective  length  scale,  u^ms/e,  is  assumed.  Therefore,  it  is  of  interest  whether  the  advective 
length  scale  is  appropriate  for  the  large  scales  of  motion  in  the  current  simulations. 
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Figure  3.G:  (a)  (11%)  vs.  Re^  (b)  v  (S2)  /  (s)  (lower  plot)  for  F\j  =  2  (•)  :  Fl  =  4  (■) 


There  is  a  strong  theoretical  argument  that  L/,  should  be  related  to  Uft  (or  urms)  and 
c.  With  appropriate  spatial  and  temporal  averaging,  advection  can  be  expected  to  balance 
the  viscous  dissipation  rate  with  an  advective  time  scale  ta  ~  L}Ju}t.  If  this  is  the  case  then 


d_\ 
dt  2 


uh 


£  ~ 


Lh 


(3.4) 


or  L/j  ~  uh/6'  While  this  analysis  is  generally  accepted  for  theoretical  estimates  of  the 
length  scale,  for  practical  application  in  simulations  or  laboratory  experiments  it  may  not 
be  sufficient  for  several  reasons.  First,  advection  of  horizontal  kinetic  energy  only  balances 
viscous  dissipation  when  the  horizontal  and  vertical  kinetic  energies  and  the  potential  energy 
are  in  equilibrium  or  when  a  carefully  chosen  interval  in  time  is  chosen  for  averaging.  Second, 
horizontal  kinetic  energy  is  convected  in  the  vertical  direction  and  u\/e  varies  by  several 
orders  of  magnitude  between  planes  of  high  and  low  shear.  As  a  consequence,  Lfx  depends 
strongly  on  the  definition  of  the  spatial  averages  used  in  its  computation. 

To  illustrate  the  difficulty  in  using  u\/e  for  analyses  of  the  current  simulations,  we 
consider  an  average  advective  length  scale 


(Lo)h  — 


(Oh 

(€)  II 


(3.5) 


computed  for  the  planes  of  maximum  shear  and  plot  it  versus  time  in  Fig.  3.7  for  each 
of  the  simulation  cases.  Here  (-)7/  denotes  a  planar  averaged  quantity.  Averaging  in  this 
manner  is  justified  by  the  fact  that  most  of  the  kinetic  energy  and  most  of  its  dissipation 
rate  is  associated  with  the  planes  of  maximum  shear,  and  it  is  in  the  vicinity  of  those  planes 
where  turbulence  occurs.  If  the  spatial  average  is  computed  over  the  entire  domain  then  the 
magnitude  of  (La}}1  is  smaller  but  the  trend  in  time  is  similar.  Note  that  with  the  the  size 
of  the  numerical  domain  and  the  fact  that  the  domain  is  periodic,  the  maximum  permissible 
value  for  the  horizontal  length  is  2n.  (La}H  is  therefore  unphysically  large  for  most  of  the 
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Figure  3.7:  Advective  length  scale  from  (3.5). 


simulation  cases.  Furthermore,  (La) varies  significantly  between  the  simulations  and  does 
not  grow  continually  in  time,  neither  of  which  is  consistent  with  information  about  the 
large  scales  of  the  flow  gained  from  the  stream  function  and  energy  spectra.  It  is  concluded 
that  ( La)jj  is  not  useful  for  estimating  the  size  of  the  horizontal  motions  in  the  current 
simulations. 

In  order  to  arrive  at  an  appropriate  length  scale  for  use  with  the  F^Re^  scaling,  we 
consider  again  the  underlying  physical  justification  for  that  scaling  (i.e.,  that  turbulence  is 
triggered  by  vertical  shearing  between  quasi-horizontal  vortices).  Based  on  this  model,  the 
appropriate  length  scale  is  the  size  of  the  horizontal  vortices.  From  the  horizontal  stream 
function  in  Fig.  3.8,  it  is  apparent  that  the  size  of  these  vortices  increases  from  one  quarter  to 
one  half  the  size  of  the  computational  domain  between  t  =  0  and  t  =  20.  A  straightforward 
approach  to  computing  Ly  is  to  relate  it  to  the  average  of  the  autocorrelations  of  u  in  the 
^-direction  and  of  v  in  the  y  direction, 

R(r)  -  1  M-r  +  r)?t(J’))/.  1  Mi/ +  rMy))/t  n6) 

2  (u)h  2  (v)h 

R(r )  is  plotted  in  Fig.  3.8  for  case  F2R32  at  two  different  times. 

The  horizontal  length  scale  is  defined  in  terms  of  the  autocorrelation  function  as 

{Lh)}f  =  r  where  R(r )  =  0  .  (3.7) 

In  Figure  3.9,  (Lfx) H  from  (3.7)  is  plotted  versus  time  for  cases  with  high  and  low  (Ref,)^. 
The  length  scale  increases  with  time  as  expected  and  never  exceeds  the  limiting  value  of 
27r.  Unlike  with  ( La)H ,  there  is  little  difference  between  the  cases,  which  is  consistent  with 
information  gained  from  the  streamfunction  and  from  spectra  that  indicate  little  difference 
in  the  evolution  of  the  large  scales  between  the  different  cases.  Furthermore,  (L}x)  H  grows 
monotonically  in  time  as  expected. 

Also  shown  in  Fig.  3.9  is  the  RMS  horizontal  velocity,  the  horizontal  Froude  number, 
and  the  horizontal  Reynolds  number.  Since  U}x  decreases  and  (Lfx) increases  with  time,  F/* 
decreases  in  time.  This  data  supports  the  theoretical  argument  by  Riley  and  de  Bruyn  Kops 
[2003]  that  stably  stratified  flows  with  no  energy  input  will  eventually  enter  the  strongly 
stratified  regime  even  if  the  initial  Froude  number  is  much  greater  than  unity.  Note  also 
from  the  figure  that  the  Reynolds  number  increases  as  the  flow  evolves,  reminiscent  of  two 
dimensional  turbulence. 
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Re=3200  F=2  t=0  Re=3200  F=2  t=20 


Figure  3.8:  Horizontal  stream  function  and  corresponding  autocorrelation  function  R(r)  for 
a  plane  of  maximum  shear  at  time  t  =  0  (left)  and  t  =  20  (right).  R(r)  is  defined  in  (3.6) 
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Figure  3.9:  (uh)^,  F and  Re/,  versus  time  for  two  cases  with  low  (Ret,)//  and  two 

cases  with  high  (Reb)//. 
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Figure  3.10:  F^Re/t  versus  (Ri)^.  The  solid  line  is  the  least-squares  linear  fit  to  the  log  of 
the  quantities.  The  circles  and  squares  represent  the  Fl  =  2  and  Fl  =  4  cases,  respectively. 


3.4.4  Parameterization  of  Turbulence 


Having  determined  a  suitable  length  scale  from  which  to  compute  F/t  and  Re/t  in  §3.4.3,  we 
now  consider  the  F^Re/,  scaling,  its  relationship  to  Ren,  and  its  utility  for  parameterizing 
turbulence  in  the  simulated  flows.  Recall  that  the  analysis  in  Riley  and  de  Bruyn  Kops 
[2003]  relates  the  Richardson  number.  Hi.  to  l/F^Re^  so  that  if  F^Re/j  >  0(1)  then  Ri  can 
be  expected  to  be  order  one  or  less.  Thus,  a  flow  with  F^Re/,  >  0(1)  will  be  susceptible  to 
Kelvin-Helmholtz  instabilities  and  turbulence  can  be  expected  to  develop.  To  investigate 
this  argument,  we  examine  the  assumption  that  Ri  1/F^Re/j  by  plotting  F^Re/,  versus 
(Ri)//  in  Fig.  3.10  for  all  of  the  simulation  cases.  Here  (Ri)//  is  the  planar  averaged  gradient 
Richardson  number: 


ZJL  (dp  +  <M\ 
p0  yaz  '  dz  J 

'WTW2' 


(3.8) 


The  results  are  encouraging  because  the  relationship  Ri  ^  1/F^Re/*  holds  over  a  two  decade 
range  of  values.  There  is  some  scatter  in  the  data,  but  there  is  no  tendency  to  deviate 
from  the  relationship  even  at  the  extreme  values  of  Ri.  In  fact,  the  stronger  relationship 
Ri  «  1/F^Re/j  is  justified  for  the  current  simulations.  The  fact  that  F^Re^  is  a  good 
estimate  for  the  Richardson  number  combined  with  the  conclusion  from  §3.4.1  that  shear 
instabilities  are  the  major  cause  of  turbulence  in  the  simulations  leads  us  to  conclude  that 
is  a  useful  parameter  for  predicting  if  turbulence  will  occur  in  these  flows. 

As  noted  in  the  introduction  of  this  paper,  the  buoyancy  Reynolds  number  can  be 
written  in  terms  of  a  Reynolds  number  and  the  square  of  a  Fronde  number.  This  leads  to  the 
question  of  whether  (Re^)^  and  FfjXo^  are  related  quantities,  and  in  Fig.  3.11  one  is  plotted 
versus  the  other  for  all  the  simulation  cases.  It  is  evident  that  the  two  parameterizations 
are  equivalent  to  within  the  scatter  of  the  data  and  an  order  one  multiplicative  constant, 
and  that  Riley  and  de  Bruyn  Kops  [2003]  arrived  at  a  new  physical  justification  for  why  Re*, 
has  proven  useful  for  parameterizing  turbulence  in  stably  stratified  flows.  This  alternative 
justification  may  help  in  understanding  the  conditions  under  which  Reb  can  be  used  to 
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Figure  3.11:  (Ret>)//  versus  F^Re/*.  The  solid  line  is  the  least-squares  linear  fit  to  the  log  of 
the  quantities.  The  circles  and  squares  represent  the  Fl  =  2  and  F^  =  4  cases,  respectively. 


parameterize  turbulence,  and  how  the  parameterization  might  be  improved. 

In  addition  to  providing  an  alternative  physical  explanation  for  the  occurrence  of  turbu¬ 
lence  under  stable  stratification,  the  F^Re^  scaling  has  several  attractive  features  compared 
with  Rob-  At  the  theoretical  level,  it  involves  two  dimensionless  groups,  which  is  the  number 
predicted  for  this  problem  by  the  Buckingham  Pi  theorem.  This  suggests  that  turbulence  pa¬ 
rameterization  be  considered  in  the  two  dimensional  Froude- Reynolds  number  space  rather 
than  in  the  one  dimensional  domain  of  a  modified  Reynolds  or  modified  Froude  number. 
At  the  practical  level  for  numerical  and  laboratory  experimentalists,  F2Re  can  be  used  a 
priori  to  estimate  if  a  flow  can  be  expected  to  be  sufficiently  turbulent  to  be  interesting  for 
understanding  oceanic  and  atmospheric  flows. 


Chapter  4 

Vortex  Street  Simulations 


This  chapter  contains  results  of  numerical  simulations  initialized  with  a  von  Karman  vortex 
street  and  a  hyperbolic  tangent  density  profile.  An  overview  of  the  simulations  are  given, 
followed  by  the  kinetic  and  potential  energy  equations  for  non-uniform  density  stratification, 
numerical  considerations,  and  simulation  results.  Emphasis  is  be  placed  on  the  effect  of 
assuming  a  density  gradient  that  is  uniform  in  height  when  it  may  not  be  uniform,  such 
as  in  a  thermohaline  staircase.  Note  that  in  keeping  with  the  convention  throughout  this 
document,  dimensional  quantities  are  denoted  with  a  tilde  (7),  nondimensiorial  quantities 
have  no  marking. 


4.1  Overview 


High  resolution  direct  numerical  simulations  (DNS’s)  of  a  perturbed  von  Karman  vortex 
street  were  performed,  simulating  the  resulting  flow  of  an  object’s  wake  (Figure  4.1).  The 
initial  conditions  consist  of  three  vortex  pairs  and  low-level  noise.  Each  simulated  flow  was 
conducted  with  no  ambient  shear  and  density  stratification  that  was  held  constant  in  time, 
to  represent  the  persistent  stratification  naturally  found  in  a  thermohaline  staircase  (§1.2.5) 
or  atmospheric  layer  transition  (§1.1.1).  The  wake  has  zero  mean  velocity,  which  is  similar 
to  that  generated  by  a  self-propelled  object  (although  Meunier  and  Spedding  [200G]  note 
that  it  is  very  difficult  to  obtain  a  truly  momentumless  wake  in  a  stratified  fluid).  Each 
vortex  was  initialized  with  the  following  velocity  profile  [de  Bruyn  Kops  et  al.,  2003]: 


V#  =  U—  exp 
r  m 


(4.1) 


where  U  is  the  initial  velocity  scale,  fm  is  a  radial  length  scale,  &u  is  the  vertical  length 
scale,  and  f  =  yjx2  T  y2  and  z  are  the  radial  and  vertical  position.  The  separation  distances 
between  vortex  centers  in  the  x  and  y  directions  were  sx  =  2 rm  and  sy  =  1.5fm. 

In  creating  the  initial  flow  condition,  noise  was  applied  to  the  vortex  horizontal  and 
vertical  length  scales  fm,  and  vortex  separation  distances  sx,  sy ;  each  was  randomly 
perturbed  up  to  5%  of  its  corresponding  scale.  For  example,  the  vertical  scale  for  each 
vortex  was  calculated  as  5u  +  0.05A<5^,  and  the  y  positions  for  the  positive  vortices  were 
determined  as  Ly/ 2  +  sy/2  +  O.OSASy,  where  A  is  a  [-1  1]  uniformly  distributed  random 
number  and  Ly  is  the  spanwise  (y)  domain  width. 
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Figure  4.1:  Center  plane  of  vortex  street  initial  condition.  Arrow  length  represents  fluid 
velocity. 


The  ambient  density  p(z)  was  imposed  with  a  hyperbolic  tangent  vertical  profile: 


p(5)  = 


from  which  the  density  stratification  dp(z)/dz  is  obtained: 


dpjz) 

dz 


(4.2) 


(4.3) 


where  A p  =  ptop  —  pbottom  is  the  difference  in  density  between  the  top  and  bottom  of  the 
numerical  domain,  5  is  the  vertical  position  and  8P  is  a  vertical  length  scale  of  the  density 
profile.  Since  the  stratification  is  a  function  of  vertical  height,  the  buoyancy  frequency  N  will 
also  be  a  function  of  height.  It  will  be  convenient  to  define  a  global  buoyancy  frequency  N  to 
describe  the  average  ambient  density  stratification  Ap/Az  =  {ptoP~  ^bottom) /{hop- bottom), 
as  well  as  a  local  buoyancy  frequency  N(z): 


Zl2 

N 


N2(*) 


a_^p 

A)  Az 

a  dji(z) 

Po  dz 


(4.4) 

(4.5) 


4.2  Theoretical  Considerations 

4.2.1  Kinetic  Energy 

The  kinetic  energy  equation  is  formed  by  taking  the  dot  product  of  the  velocity  vector  v 
and  the  momentum  equation  (2.46b)  (in  which  the  Boussinesq  approximation  is  made  and 
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gravity  is  the  only  body  force).  This  yields  (in  indicial  notation): 


dEk  .  dEk  _  . 

dt  U 1  dij  Ul  dxi  po 


dp  9  ~~  ,  d(uiTji)  ~ 

'  =  Pt‘3  +  T^d£j 


dxi 


(4.6) 


Here  Ek  =  1/2  (u«Ut)  is  the  kinetic  energy  per  unit  mass,  aq  is  the  direction  vector,  and 
Tji  is  the  viscous  stress  tensor  for  a  Newtonian  fluid.  Equation  (4.6)  can  be  broken  into 
horizontal  (i  =  1,2)  and  vertical  (i  =  3)  contributions  (denoted  with  subscripts  II  and  V 
respectively): 


dEk 

dt 

dEn 

dt 

dEv 

dt 


dEp  dEy 


dt  '  dt 

=  -Til  —  Ph  +  Wjj  —  £h  -  4> 

=  -Tv  -  Pv  —  B  +  W\.  -  ev  +  4> , 


(4.7) 

(4.8) 

(4.9) 


where 


Eh=^ 

Tu  =  v  •  VEh 
Ph  =  uh  •  V  HP 

m  r  3(uh  Tjh) 

WH  =  ()t.j 

Eli  =  20c.nje.iij 

<p  =  2t>(ei:jf  13  +  623**23) 


Ev=  ^ 

Ty  =  v  •  VEv 

Pv  =  uv§ 

W\  = 

Ey  =  2 OayjCyj 

B  =  H-pP 


In  the  above  equations  v  —  (u,v7w)  is  the  velocity  vector.  Uh  =  (u.u)  is  the  horizontal 
component  of  the  velocity  vector,  e  =  2 ueijCij  is  the  dissipation  rate  of  kinetic  energy. 
0  =  2i/(ci3ri3  +  ^23^23)  is  a  coupling  term  between  vertical  and  horizontal  kinetic  energy, 


ea 


and  are  the  symmetric  and  anti-symmetric  tensors  of 


diii 

()x-i 


eij  ~ 


Vij  — 


Nondimensionalization 


nondirnensional  form  as: 


dE„ 

dt 

dEv 

dt 


diii  diij  \ 

;  Ci  .z,  ) 

(4.10) 

t'Xj  OXi  J 

diu  du  7  \ 

(4.11) 

*  j 

dij  dii  J 

equations  (4.8)  and  (4.9) 

can  be  written  in 

/  +  Wh  -  eh  ~<t> 

(4.12) 

r  —  B  +  1+ \/  — 

(4.13) 

where 
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Eh  = 

Tn  =  v  ■  VEh 


PlI  =  Uh  •  V HP 

W '*  = 

£n  =  ikeiUeHj 

0  =  ^:(ei3?'l3  +  e237’23) 


Ey  = 

Ty  =  V  •  V  Ey 

Pv  =  UV'jl 

£v  =  m;evjGvj 

B  =  (R)V. , 


and  Rer,  Fr,  Sc  are  defined  in  (4.28).  Further  discussion  of  each  term  is  given  in  the  results 
section  of  this  chapter. 


4.2.2  Available  Potential  Energy  for  Non-uniform  Density  Stratification 

Potential  energy  in  a  geophysical  setting  usually  involves  the  concepts  of  available  and 
background  potential  energy,  first  suggested  by  Lorenz  [1955].  He  noted  that  in  order  to 
convert  the  total  potential  energy  in  the  Earth’s  atmosphere  to  kinetic  energy,  the  temper¬ 
ature  needed  to  reach  absolute  zero  and  all  mass  needed  to  be  located  at  sea  level.  Such 
conditions  cannot  readily  occur.  (It  is  estimated  that  potential  energy  makes  up  25%  of 
the  total  energy  (internal  +  potential  +  kinetic)  in  the  Earth’s  atmosphere,  while  only  2% 
is  kinetic  energy  [Gill,  1982,  pg.81]).  Instead,  the  potential  energy  Ep  that  is  available  for 
conversion  to  kinetic  energy  is  said  to  be  the  result  of  any  deviation  from  a  background 
(or  rest)  potential  energy  E b.  E b  is  a  state  that  would  exist  if  the  fluid  was  adiabatically 
redistributed  (i.e.,  no  heat  transfer)  to  a  minimum  energy  state.  The  available  potential 
energy  is  the  total  potential  energy,  P,  minus  the  background  potential  energy: 

Ep  =  V-Eb .  (4.14) 

Initially,  this  might  seem  like  a  straightforward  method  for  obtaining  the  available  po¬ 
tential  energy  in  a  system.  However,  while  V  is  typically  defined  as  Jf  J  ptgz  dx  dy  dz ,  (where 
Pt  —  Po  +  p(^)+P  is  the  total  density,  equal  to  the  sum  of  reference,  ambient  and  fluctuating 
components),  several  methods  of  obtaining  Ep  exist.  Typically,  adiabatic  redistribution  is 
performed  by  sorting  the  density  field  so  that  the  highest  density  parcels  are  in  the  lowest 
vertical  position.  The  redistributed  density  field,  p*,  is  then  used  to  calculate  background 
potential  energy  Eb  =  fff  p*gz  dx  dy  dz  [e.g.,  Staquet,  2000,  Winters  et  al.,  1995].  Another 
related  method  to  calculate  Eb  is  by  “Thorpe  reordering,”  [Thorpe,  1977],  in  which  the 
density  is  also  adiabatically  redistributed  via  sorting,  but  in  this  case  the  distance  (absolute 
or  rms)  each  fluid  element  traversed,  f/,  to  obtain  the  minimum  energy  state  is  used  to  cal¬ 
culate  the  background  potential  energy  Jff  ptg(z  —  y)dxdydz  [Smyth  and  Mourn,  2000a, 
Smyth  et  al.,  2001].  A  third  method  of  obtaining  £),,  suggested  by  Tseng  and  Fcrziger 
[2001],  involves  taking  the  probability  density  function  (PDF)  of  pt)  which  can  be  thought 
of  as  a  method  of  sorting  the  field  into  a  minimum  energy  state.  With  the  PDF  of  p*,  the 
vertical  position  of  each  parcel  in  a  minimum  energy  state  can  be  found,  then  integrated 
over  the  domain  height  to  find  Eb. 

While  the  above  methods  for  obtaining  Eb  rely  on  sorting  methods,  thus  avoiding  the 
need  for  derivatives  of  density  stratification,  there  are  several  drawbacks.  First,  as  pointed 
out  in  Winters  et  al.  [1995],  sorting  methods  can  only  provide  an  estimate  of  background 
potential  energy,  since  the  vertical  position  is  discretized  by  the  numerics,  where  in  a  physical 
system  no  such  discretization  exists.  Second,  since  simulations  arc  now  becoming  large  (it 
is  not  uncommon  for  simulations  to  be  the  order  of  1  billion  gridpoints),  sorting  or  creating 
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a  PDF  of  a  large  field  can  take  too  long  to  be  practical  if  values  are  to  be  obtained  while 
the  simulation  is  running.  Lastly,  to  gain  insight  into  Eulerian  energetics,  it  is  desirable  to 
calculate  Ep  on  a  local  basis  (i.e.,  at  each  grid  point),  whereas  the  above  methods  are  only 
defined  on  a  volume  average  basis.  A  method  for  computing  Ep  on  a  local  basis  is  derived 
below. 

Available  potential  energy  for  uniformly  stratified  flows  is  usually  defined  locally  as  [e.g., 
Gill,  1982,  p.140]: 


EP  = 


1  A  ~2 

2  dpjdiz 


(4.15) 


where  p(x,  y,z,  t)  is  the  fluctuation  from  the  undisturbed  density  p,  and  dp/dz  is  constant 
in  space  and  time.  When  density  stratification  is  non-uniform,  the  definition  of  potential 
energy  is  more  complicated  than  for  uniform  stratification,  and  (4.15)  will  not  be  accurate. 
This  complication  arises  from  the  fact  that  the  derivative  of  the  stratification  is  non-zero 
(i.e..  d2p(z)/dz 2  ^  0)  [Holliday  and  McIntyre,  1981] .  A  method  to  obtain  an  accurate 
expression  for  potential  energy  in  an  incompressible  fluid  with  a  non-uniform  density  strat¬ 
ification  was  put  forth  by  Holliday  and  McIntyre  [1981],  who  begin  by  defining  available 
potential  energy  as  the  integral  of  the  displacement  of  a  fluid  particle  from  its  undisturbed 
state  (0, 

Ep(z,0  =  ~  j  ~C)dC.  (4-16) 

where  (:)  denotes  a  dummy  integration  variable.  Since  the  numerical  simulations  involve 
a  field  of  density  fluctuations  p(x,  y,  2.  t).  it  is  advantageous  to  write  Ep  in  terms  of  p(z) 
and  p(z,  ?y,  z,  £),  rather  than  z  and  Conversion  to  p(z)  and  p(.r,  j),  z ,  t)  is  possible  because 
p(:r,  y .  z,  t )  contains  Lagrangian  information  required  to  compute  (4.1  G^.  Converting  to  p(z) 
and  p(x,y,zj)  can  be  achieved  by  first  defining  a  potential  function  4>{}  as 


${^(^)}  =  gz. 


(4.17) 


Provided  the  undisturbed  stratification  p(z)  is  stable  everywhere,  (4.1G)  can  be  written  in 
the  following  form  [Holliday  and  McIntyre,  1981  : 


=-/ 


4>{p(z)  +  p]  -  4>{p{z)}  dp, 


(4.18) 


where  p  is  a  dummy  integration  variable  and  the  spatial  dependence  (x,  i),  z)  of  p  is  implied. 
For  the  specific  hyperbolic  tangent  p{z)  given  in  (4.2),  (4.17)  can  be  shown  to  be 


4>]p(z)}  —  —  gSp  arctanh 


(4.19) 


where  p  =  Ap/2,  obtained  from  (4.2),  is  introduced  simply  for  notational  purposes.  Substi¬ 
tuting  (4.19)  into  (4.18)  yields  the  following  expression  for  Ep  (where  the  spatial  dependence 
of  p,  p  is  implied): 


Ep  =  gSp 


(p  +  p)  arctanh 


ep 


e  -  pp-  p 


tP  l+  f  111 


Q  -  (P  +  p) 


2  \  I 


Q2  ~  p2 


The  time  derivative  of  Ep  can  be  found  by  use  of  the  chain  rule: 


0EP  _  0Ep  dp 
~df  =  ~Wdt' 


(4.20) 


(4.21) 
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Note  the  dEp/dp  is  the  opposite  of  the  integral  (4.18).  Thus,  for  the  hyperbolic  tangent 
stratification  studied  here,  dEp/di  is 

Op 


=  p<Larctanh 


QP 


dt  p  \g2  —  pp  —  p2  J  dt 

Substitution  of  the  density  trasoport  equation  (2.46c)  into  (4.22)  yields: 


OE, 


dt 


J-  =  <j<Larctanh 


QP 


Q2  -pp-  p2 


-V  -Vp-w^l  +  DV2p)  , 
dz 


(4.22) 


(4.23) 


where  the  material  derivative  in  (2.46c)  has  been  expanded,  and  the  symbol  D  replaces  k 
for  the  mass  diffusivity. 

Equation  (4.23)  can  be  written  in  shortened  notation,  similar  to  the  kinetic  energy 
equations  (4.8)  and  (4.9): 


OE, 


£  =  -Tp  +  B-x 


where 


Tp  —  gfip 
B  =  -gSp 
x  =  -gtpD 


dt 

arctanh 

arctanh 


(4.24) 


gp 


Q2  ~  PP  ~  P2 


QP 


Q2  -  pp-  p2 


QP 


v  •  Vp 

-  dp{z) 
dz 
2  - 


arctanh  — - — - v  p 

'  g2  —  pp  —  p2  ' 


An  interesting  question  that  arises  is  how  does  (4.20)  differ  from  the  potential  energy 
defined  by  (4.15).  This  can  be  seen  by  expanding  (4.20)  in  a  Taylor  series  of  about  p  =  0. 
With  some  algebraic  manipulation,  the  first  3  terms  of  the  Taylor  series  can  be  shown  to 
be: 


Ep  = 


1  a 


l 


-2 


1  9 


p{z) 


-3 


l  (J  Q  +  3  p(z) 


2  po  dp(z)/dz^  3  po  Spg(dp(z)/dz)2^  12  po  gSp(dp(z)/dz)3^ 

In  the  limit  of  linearized  motion  (which  is  typically  how  (4.15)  is  derived),  terms  of  order 
greater  than  p2  are  assumed  small  and  neglected,  reducing  the  above  equation  to  (4.15). 


Nondimensionalization  Using  the  scaling  outlined  in  §4.4.1  below,  (4.24)  can  be  written 
in  nondiineiisional  form  as: 

°EP  - Tp  +  B-x ,  (4.26) 


where 


2n\  J 


ef  -  I  Tr)  sp 

*  ■  l'g)> 

B  =  (S)  <5 


X 


dt 


( p(z )  +  p)arctanh 


QP 


arctanh 


QP 


Q2  ~  PP-  P2 


arctanh 


QP 


Q2-  PP-  P2 


=  (^)2  — 
\FrJ  RerSc 


arctanh 


,  e , 

9  -  -2  /  “*"  O  11 

r'  ~  pp  —  p  !  2 


v  •  Vp 

dpjz ) 

dz 

V2p 


g2  -  (p(j 0  -  p) 

p2  —  p2 


2  \  1 


W- 


QP 


Q2  -PP-  P2 


and  Rer,  Fr,  Sc  are  defined  in  (4.28)  below. 
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4.3  Numerical  Considerations 

Naturally,  it  is  not  sufficient  to  perform  numerical  simulations  without  first  performing 
basic  numerical  checks  for  adequate  spatial  resolution,  temporal  resolution,  and  domain  size. 
Spatial  resolution  includes  both  large  and  small  scales.  If  the  simulation  is  not  spatially 
resolved,  energy  will  not  transfer  to  proper  scales,  usually  resulting  in  an  under-calculation  of 
kinetic  energy  dissipation  rate.  Since  the  simulations  use  a  spectral  code,  spatial  resolution 
can  be  verified  via  an  energy  spectrum.  If  the  simulation  is  not  resolved  temporally,  each 
timestep  will  yield  an  inaccurate  result,  and  this  will  affect  the  energy  balance.  As  the 
simulations  will  be  carried  out  with  a  variable  timestepping  technique,  temporal  resolution 
can  be  checked  by  integrating  the  left  hand  side  of  the  kinetic  and  potential  energy  equations 
(4.8),  (4.23).  Finally,  the  vertical  domain  size  may  have  an  effect  on  the  simulation.  If  the 
domain  is  too  small,  the  boundaries  may  adversely  affect  the  simulation.  Each  of  these 
numerical  considerations  is  examined  below. 


4.3.1  Small  Scale  Resolution 

Small  scale  motions  must  be  resolved  in  a  DNS  to  properly  simulate  the  dissipation  rate  of 
kinetic  and  potential  energy.  One  method  of  verifying  the  resolution  of  the  flow  is  by  means 
of  the  kinetic  energy  dissipation  rate  (e)  spectrum.  Since  €  is  a  small  scale  process,  a  well 
resolved  simulation  will  show  the  magnitude  of  the  kinetic  energy  dissipation  rate  spectrum 
to  increase  with  wavenumber  (decreasing  length  scale),  indicating  e  is  occurring  at  small 
scales.  Since  the  How  is  initialized  with  maximum  energy  on  the  center  horizontal  plane, 
this  would  be  a  good  place  to  check  the  simulation  resolution,  as  it  is  a  likely  place  for  the 
simulation  to  become  under-resolved.  Figure  4.2  contains  a  plot  of  the  (nondimensional) 
kinetic  energy  dissipation  rate  spectrum  of  the  center  plane  at  t  =  10,  with  Fr  =  2.75. 
Rer  =  19200,  £  =  0.01,  which  is  representative  of  all  vortex  street  simulations  performed 
in  this  study.  In  Figure  4.2,  e  increases  with  wavenumber  up  to  kx  =  112,  where  the  anti¬ 
aliasing  filter  causes  a  sudden  drop  in  magnitude.  This  result  suggests  that  the  simulation 
is  well-resolved.  Several  other  horizontal  planes  (including  planes  of  maximum  shear  seen 
in  Figure  4.11)  were  examined  with  similar  results. 

4.3.2  Temporal  Resolution 

Within  the  Boussinesq  approximation,  the  mechanical  energy  equation  and  momentum 
equations  are  not  independent.  Since  the  momentum  equations  are  numerically  integrated, 
the  energy  equations  can  provide  a  check  of  the  temporal  resolution.  In  order  to  verify 
temporal  resolution,  the  left  hand  side  of  the  volume  averaged  kinetic  and  and  potential 
energy  equations  (eqs.  (4.6),  (4.23))  are  integrated  in  time  using  the  trapezoid  rule  and 
compared  to  the  directly  computed  energy  values.  Figure  4.3  (a)  contains  plots  of  the 
integrated  kinetic  energy  equation  vs.  calculated  energy  for  simulation  £  =  0.01,  Fr  =  2.75, 
Rer  =  19200,  while  Figure  4.3(b)  contains  a  plot  of  the  relative  error  between  the  two. 
The  same  quantities  for  potential  energy  are  shown  in  Figures  4.4  (a),  (b).  In  each  plot, 
there  is  good  agreement  between  the  calculated  and  integrated  energy  values.  Early  in 
each  simulation,  the  relative  error  is  «  10-3.  The  relative  error  is  seen  to  increase  in  time, 
and  can  be  explained  by  noting  that  as  each  term  in  integrated,  the  cumulative  error  will 
increase,  which  will  gradually  increase  the  relative  error.  In  addition,  data  was  available 
in  ^  0.25  nondimensional  time  units,  which  is  a  large  timestep  to  be  integrating  over. 
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Figure  4.2:  (nondimensional)  x  direction  kinetic  energy  dissipation  rate  spectrum  (sn)  for 
the  center  plane  of  £  =  0.01,  Fr  =  2.75,  Rer  =  19200. 


Also,  the  integration  was  performed  with  a  first  order  trapezoid  method  for  convenience,  as 
opposed  to  a  more  accurate  alteration  method  like  3rd  order  Adams  Bashforth.  When  the 
simulation  was  run  for  7  consecutive  tiinesteps,  the  relative  error  between  the  integrated 
and  actual  energy  values  is  «  10  *  (Figure  4.5).  Again,  the  increase  in  error  is  due  to 
the  cumulative  error  added  for  each  tiinestep  integration.  These  results  suggest  that  the 
simulation  is  temporally  resolved. 

4.3.3  Vertical  Domain  Size 

When  analyzing  stratified  flows,  the  behavior  of  vertical  motions  is  a  topic  of  interest.  As 
such,  it  is  important  to  consider  the  effect  of  the  vertical  domain  size  on  the  simulated 
flows.  Since  the  simulations  are  run  with  periodic  boundary  conditions,  it  is  similar  to 
having  an  identical  flow  domain  at  each  boundary.  Thus,  flows  that  occur  at  one  edge  of 
the  domain  can  affect  the  flows  at  the  opposite  boundary.  The  effect  of  this  interaction  can 
be  seen  in  Figure  4.6,  which  contains  plots  of  the  horizontal  (xy)  planar  averaged  horizontal 
and  vertical  kinetic  energies  (FT/,  Ey)  at  t  =  10  for  simulations  with  nondimensional 
computational  domain  height  Lz  =  Lz/fm  =  ±1.5  and  ±3.  Note  that  Eh  is  similar  for  each 
simulation  except  for  the  vertical  boundaries,  where  the  energy  increases  at  the  boundary  for 
Lz  =  ±1.5  as  opposed  to  when  Lz  =  ±3.  In  the  vertical  kinetic  energy  plot,  a  large  double 
peak  exists  above  and  below  the  centerline  with  Lz  =  ±1.5,  but  decreases  significantly  in 
magnitude  with  Lz  =  ±3.  Also,  an  internal  wave  can  be  seen  in  Ey  (and  to  a  lesser  extent 
in  Eh)  that  extends  to  the  boundary  with  Lz  —  ±3,  but  does  not  exist  with  the  smaller 
domain  height.  Support  for  internal  wave  formation  can  be  seen  in  Figure  4.18  which  is 
consistent  with  the  fact  that  internal  waves  can  not  exist  when  there  is  no  stratification. 
The  differences  in  Lz  can  be  attributed  to  the  periodic  boundary  condition,  where  for  the 
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Figure  4.3:  Kinetic  energy  equation  balance  (left)  and  relative  error  between  integrated  and 
directly  computed  kinetic  energy  (right).  (•)  represents  the  value  of  energy  at  a  particular 
time,  (-)  is  the  time  integration  of  the  right  hand  side  of  (4.G).  Here  F,  =  2.75.  Rer  =  19200, 
£  =  0.01 


Figure  4.4:  Potential  energy  equation  balance  (left)  and  relative  error  between  integrated 
and  directly  computed  potential  energy  (right).  (•)  represent  the  value  of  energy  at  a 
particular  time,  (-)  is  the  time  integration  of  the  right  hand  side  of  (4.23).  Here  Fr  =  2.75, 
Rer  =  19200,  £  =  0.01. 
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Figure  4.5:  Relative  error  between  integrated  and  computed  kinetie  and  potential  energies 
for  7  consecutive  timesteps,  simulation  Fr  =  2.75,  Rer  =  19200,  £  =  0.01. 


smaller  Lz  the  wave  interacts  more  strongly  at  the  boundaries  than  the  larger  Lz.  For  all 
simulations  Lz  was  chosen  to  be  ±3  to  minimize  the  interaction  between  boundaries. 


4.4  Methodology 

4.4.1  Equations  of  Motion 

The  numerical  method  used  in  the  vortex  street  simulations  is  the  same  as  the  Taylor-Green 
simulations  in  Chapter  3.  Speeifieally,  the  fields  are  assumed  to  satisfy  the  Navier-Stokes 
equations  subject  to  the  Boussinesq  approximation  (§2.1.5),  and  are  solved  using  a  pseudo- 
speetral  technique  with  periodic  boundary  conditions  and  the  pressure-projeetion  method. 
Taking  U  as  the  velocity  scale,  frn  as  a  length  scale,  r7n\Ap/ Az\  as  a  density  scale  (where 
|  *  |  denotes  absolute  value),  fm/U  as  a  time  scale,  and  p^U2  as  a  pressure  scale  (where  po 
is  the  reference  density  value),  the  nondimensional  governing  equations  in  a  non- rotating 
frame  of  reference  are: 

V  •  v  =  0  (4.27a) 

§F  +  V'VV  =  -(??)  ^-Vp+^-V\  (4.27b) 


dp 

dt 


+  V  •  V/9  —  W 


dp(z) 

dz 


1 


RerSc 


V2p 


(4.27e) 


where  v  =  (it,  u,  w)  is  the  velocity  vector,  and  p  and  p  are  the  (nondimensional)  density 
and  pressure  deviations  from  their  ambient  values,  and  ez  is  a  unit  veetor  in  the  vertical 
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Figure  4.G:  Comparison  of  horizonal  planar  averaged  horizontal  (left)  and  vertical  (right) 
kinetic  energy  with  Lz  —  ±1.5  and  ±3.  t  =  10.  Here  Fr  =  2.75,  Re,  =  19200,  £  =  0.01. 


direction.  The  Reynolds,  Froude,  and  Schmidt  numbers  are  defined  as: 


Rer 


v 


—  27 tU 

r  r  —  ~ 

t  rn 


(4.28) 


— 2 

where  N  =  —g/po(Ap/Az)  is  the  average  buoyancy  (or  Brunt- Vaisala)  frequency,  and  V  is 
an  “effective”  mass  diffusivity,  representing  the  combined  effects  of  thermal  diffusivity  and 
either  salt  diffusivity  (ocean)  or  water  vapor  diffusivity  (atmosphere). 

In  all  vortex  street  simulations,  Fr  =  2.75  and  Rer  =  19200.  In  addition,  Sc  was  set  to 
1,  close  to  the  ratio  of  momentum  to  heat  diffusivity  in  air  (Sca*r  =  0.7),  but  far  from  the 
oceanic  Schmidt  numbers  for  temperature  (Sct  =  7)  and  salinity  (Scs  =  700).  Also,  note 
the  only  difference  between  (4.27)  and  (3.2)  is  the  allowance  of  the  density  stratification 
dp(z)/dz  to  be  a  function  of  vertical  height,  rather  than  a  constant. 


4.4.2  Momentum  vs.  Density  Vertical  Scales 

The  vertical  profile  of  both  the  velocity  (4.1)  and  density  stratification  (4.3)  are  sech2,  with 
different  vertical  scales  Sy  and  6P.  The  parameter  £  is  now  defined  which  describes  the  ratio 
of  the  wake  to  density  vertical  length  scales: 


(4.29) 


In  each  simulation,  5y  remains  fixed,  while  Sp  is  altered.  Simulations  were  performed  with 
£  ranging  from  0.01  (nearly  uniform  density  stratification)  to  4  (sharp  density  step).  The 
vertical  density  stratification  profiles  with  several  £’s  are  shown  in  Figure  4.7.  Note  that 
while  locally  dp{z)/dz  differs  for  each  £,  the  average  (nondimensional)  change  in  density 
with  height  Ap/Az  =  —  1  is  the  same. 
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Figure  4.7:  (a)  Ambient  (nondimeiisional)  density  profile  and  (b)  stratification  profile  for 
each  Note  Ap/Az  =  —  1  for  each  blit  locally  dp(z)/dz  is  different. 


Notation 

€ 

dp 

Su 

Fr 

Iler 

Lx 

Ly 

Lz 

Nx 

Ny 

Nz 

F2.75R24 

0.01 

G9.4 

0.694 

2.75 

2400 

12 

6 

6 

512 

256 

256 

F2.75R48 

0.01 

69.4 

0.694 

2.75 

4800 

12 

6 

6 

512 

256 

256 

F2.75R192 

0.01 

69.4 

0.694 

2.75 

19200 

12 

6 

6 

1024 

512 

512 

£  =  0.01 

Same  as  F2.75R192 

£  =  0.5 

0.5 

1.39 

0.694 

2.75 

19200 

12 

6 

6 

1024 

512 

512 

*  =  1 

1 

0.694 

0.694 

2.75 

19200 

12 

6 

6 

1024 

512 

512 

S  =  2 

2 

0.347 

0.694 

2.75 

19200 

12 

6 

6 

1024 

512 

512 

£  —  4 

4 

0.174 

0.694 

2.75 

19200 

12 

6 

6 

1024 

512 

512 

No  Strat 

N/A 

N/A 

0.694 

oo 

19200 

12 

6 

6 

1024 

512 

512 

Table  4.1:  List  of  vortex  street  simulations 


4.5  Simulation  Results 

In  this  section  results  are  presented  for  vortex  street  simulations.  Comparisons  will  be  made 
with  different  Rer  for  £  =  0.01  (yielding  a  near  uniform  density  stratification,  see  4.7)  a s 
well  as  for  all  £  simulations,  with  Fr  —  2.75  and  Rer  —  19200.  Table  4.1  contains  a  list  of  the 
vortex  simulations  performed,  including  Froude  and  Reynolds  number  as  defined  in  (4.28), 
domain  size,  and  number  of  grid  points  NXJ  iVy,  Nz  in  the  x,  ?/,  and  z  directions.  General 
flow  characteristics  are  first  presented  followed  by  flow  energetics  and  finally  quantities  of 
interest  such  as  buoyancy  Reynolds  number  and  mixing  efficiency. 

4.5.1  General  Flow  Characteristics 

In  order  to  gain  insight  into  the  general  flow  dynamics,  a  time-series  plot  of  the  center  plane 
streamfunction  ip  as  defined  in  Riley  and  de  Bruyn  Kops  [2003]  is  shown  in  Figure  4.8  for 
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t=0  t=5 


Figure  4.8:  Contour  plot  of  center  plane  stream  function  ip  for  F2.75R192,  £=0.01.  Light 
colors  represent  positive  values  of  dark  colors  represent  negative  values. 


the  simulation  £  =  0.01,  where  light  colors  represent  positive  vortices,  dark  colors  represent 
negative  vortices.  Cases  having  other  £  values  demonstrate  qualitatively  similar  results. 
The  von  Karman  vortex  street  is  clear  at  t  =  0.  As  the  simulated  flow  evolves,  the  vortices 
interact  with  each  other,  and  vortex  pairing  occurs  by  t  =  15. 

As  the  flow  evolves,  vertical  velocity  (w)  is  generated  due  to  internal  waves  and  turbu¬ 
lence  formation.  Since  the  flow  was  initialized  with  no  vertical  velocity,  w  can  be  used  as 
an  indicator  of  turbulence  generation.  Figure  4.9  contains  plots  of  the  center  plane  ( z  =  0) 
vertical  velocity  at  several  different  times  for  simulation  F2.75R192.  At  t  —  5  a  “herring 
bone”  pattern  of  vertical  velocity  forms  where  there  is  large  shear  between  vortices.  The 
vertical  velocity  increases  in  time,  reaching  a  at  maximum  approximately  t  =  10  before 
decaying. 

In  density  stratified  flows  it  is  well  known  that  turbulence  forms  in  intermittent,  localized 
patches.  Figure  4.10  contains  plots  of  center  plane  vertical  velocity  at  t  =  10  for  simulations 
£  =  0.01  (left  plot)  and  No  Strat  (right  plot).  Simulation  £  =  0.01  is  representative  of  all 
stratified  simulations,  where  intermittent  turbulent  patches  form,  which  is  consistent  with 
prior  studies  regarding  density  stratified  flows.  In  contrast,  vertical  velocity  is  shown  to 
be  distributed  throughout  the  plane  for  the  No  Strat  simulation.  This  results  suggests  the 
stratified  simulations  are  representative  of  flows  observed  in  laboratory  and  natural  settings. 


Shear  and  Richardson  Number 

In  stratified  flows  dominated  by  vortical  modes  (such  as  those  simulated  in  this  study), 
it  is  postulated  that  horizontal  layer  decoupling  occurs,  and  the  flow  will  be  susceptible 
to  Kelvin-Helmholtz  shear  instabilities  [Lilly,  1983].  The  square  of  the  vertical  shear  of 
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Figure  4.9:  Center  plane  ( z  =  0)  vertical  velocity  for  the  same  times  as  in  Figure  4.8.  All 
figures  have  the  same  colormap  scaling.  Dark  colors  represent  negative  (downward)  velocity, 
light  colors  represents  positive  velocity. 


Figure  4.10:  Center  plane  (z  =  0)  vertical  velocity  for  £  =  0.01  (left  plot)  and  no  stratifica¬ 
tion  (right  plot).  Dark  colors  represent  negative  (downward)  velocity,  light  colors  represent 
positive  velocity. 
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horizontal  motions  is  defined  as: 


(4.30) 


Figure  4.11  contains  time  evolution  plots  of  the  horizontally  averaged  square  of  vertical  shear 
(S2) H  vs.  vertical  position  for  several  £  and  non-stratified  simulations,  where  (-)H  denotes 
horizontal  (xy)  planar  averaging.  Initially,  the  shear  profile  demonstrates  a  “bimodal” 
pattern  with  maximum  values  at  2  ±0.5.  a  result  of  the  initial  sech2(2)  velocity  profile  (note 
that  the  initial  vertical  shear  magnitude  is  very  small).  The  bimodal  shear  pattern  persists 
in  time  for  all  the  stratified  cases,  which  can  be  attributed  to  decoupling  of  horizontal  layers 
in  stratified  flows  suggested  by  Lilly  [1983].  In  contrast,  the  bimodal  pattern  disappears 
by  t  =  5  for  the  non-stratified  case,  where  decoupling  does  not  occur.  It  is  interesting  to 
note  that  the  magnitude  of  (S2)n  increases  as  £  is  increased.  This  is  due  to  the  increase 
in  local  density  stratification  with  (Figure  4.7),  which  in  turn  causes  stronger  horizontal 
layer  decoupling  and  an  increase  in  the  vertical  shear  between  layers.  Also,  when  £  =  4 
the  maximum  {S2)H  occurs  much  earlier  in  time  than  for  other  values  of  £.  This  is  caused 
by  a  rapid  increase  in  shear  at  the  intersection  of  the  unstratified  flow  regime  outside 
the  density  stratification  layer  with  the  strongly  stratified  regime.  That  is,  for  £  =  4  a 
significant  amount  of  (S2)fJ  occurs  where  there  is  no  density  stratification  (|2|  >  0.75), 
resulting  in  the  creation  of  vertical  overturning  in  this  area,  seen  in  the  planar  averaged 
vertical  kinetic  energy  ( Ey) H  (Figure  4.18(e)  below).  The  density  stratification  then  acts 
as  a  boundary  between  the  overturning  flow  outside  the  stratification  and  the  horizontal 
How  inside,  resulting  in  a  larger  amount  of  shear  to  be  generated. 

A  quantity  related  to  vertical  shear  of  horizontal  motions  is  the  gradient  Richardson 
number  Ri,  defined  as  the  ratio  of  buoyancy  to  shear  forces.  The  volume  and  planar 
averaged  Ri  are  defined  here  as: 


<Ri> 


<Ri>H 


vfJ  (s2) 

27tV2  dp(z)/dz 

Tr)  {S% 


(4.31) 

(4.32) 


where  (•)  denotes  volume  averaged  quantity.  (Ri)  is  shown  in  Figure  4.12(a).  Initially 
(Ri)  is  large,  indicating  that  there  is  insufficient  shear  to  cause  instabilities  and  turbulence. 
Then,  as  the  flow  progresses  in  time,  (Ri)  rapidly  decreases  due  to  vertical  decoupling  of 
the  horizontal  motions.  (Ri)  decreases  to  about  0.25,  supporting  the  notion  that  Kelvin- 
Helmholtz  instabilities  are  a  mechanism  by  which  turbulence  occurs.  Further  examination 
of  the  planar  averaged  Richardson  number  (Ri)//  (Figure  4.12(b))  shows  that  for  all  values 
of  £,  (Ri)//  is  significantly  less  than  0.25  for  planes  near  the  center  of  the  wake,  suggesting 
it  is  here  where  Kelvin-Helmholtz  instabilities  are  begin  generated. 

From  Figure  4.12,  it  is  apparent  that  the  evolution  of  (Ri)  and  the  profiles  of  (Ri)// 
depend  strongly  on  £.  When  £  <  1,  all  but  the  tails  of  the  wake  are  strongly  stratified  and 
the  stratification  in  the  center  of  the  wake  is  stronger  for  higher  values  of  £.  The  flows  with 
£  <  1  behave  as  expected  based  on  similar  simulations  in  uniformly  stratified  turbulence 
Riley  and  de  Bruyn  Kops  [2003].  In  particular,  the  development  of  the  flow  is  qualitatively 
independent  of  the  strength  of  the  stratification  in  the  core  of  the  wake,  but  the  time  for 
decoupling  of  the  horizontal  motions  to  occur  and  (Ri)  to  reach  its  minimum  value  increases 
slightly  with  increased  stratification.  When  £  >  1,  however,  the  wake  profile  is  wider  than 
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Figure  4.11:  Time  evolution  of  vertical  shear  of  horizontal  motion  (S2)n  for  all  £  and  no 
stratification.  Note  the  difference  in  horizontal  scale  for  plot  (e).  dp(z)/dz  is  shown  as 
bold  dashed  line,  and  is  scaled  differently  in  each  figure  for  comparison  to  (S'2)//.  Bold 
horizontal  lines  mark  where  dp(z)/dz  =  —0.01.  Since  the  entire  vertical  domain  satisfies 
this  criterion  for  plots  (a)-(b),  no  bold  horizontal  lines  are  shown. 
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Figure  4.12:  (a)  (Ri)  vs.  time,  (b)  Vertical  profile  of  (Ri)/7  vs.  height  at  t  —  10.  Note 
that  all  simulations  begin  with  (Ri)  =  25.  Also  note  the  profiles  of  Ri  with  £  >  2  are 
fundamentally  different  from  those  with  £  <  2. 


the  stratification  profile  and  the  regions  initialized  with  the  highest  shear  are  not  strongly 
constrained  in  the  vertical  direction  by  gravity,  resulting  in  earlier  less  inhibited  formation 
of  shear. 

For  the  simulations  with  £  <  1,  the  minimum  (Ri)//  is  located  near  z  —  0,  as  all  the 
wake  is  inside  the  density  stratification  layer.  In  contrast,  when  £  >  1,  some  of  the  wake  is 
now  outside  the  stratification  layer,  and  minimum  (Ri)//  occurs  at  \z\  >  1.5.  This  behavior 
of  Ri  helps  explain  the  large  vertical  kinetic  energy  shown  in  Figure  4.18  below,  particularly 
for  £  —  4. 


Horizontal  Length  Scale 


Flows  subject  to  stable  density  stratification  often  demonstrate  an  increasing  horizontal 
length  scale.  This  quality  is  evident  in  the  evolution  of  the  center  plane  streamfunction 
(0),  shown  in  Figure  4.8.  Theoretical  arguments  suggest  that  a  horizontal  length  based  on 
an  advective  time  scale  I n  ~  u^/e  can  be  calculated.  However,  in  flows  that  demonstrate 
intermittent  turbulence,  use  of  such  an  advective  length  scale  can  result  in  unrealistic  values 
of  £//  (see  §3.4.3).  Instead,  a  horizontal  length  scale  is  defined  here  based  on  the  cross¬ 
correlation  of  velocity: 


Ryx{r)  — 


{ v(x  +  r)v(x)) 


ii 


<»2> 


H 


(4.33) 


where  the  Ryx  denotes  the  correlation  of  y  direction  velocity  (v)  in  the  x  direction.  The 
horizontal  length  scale  is  then  defined  as  the  distance  r  when  Ryx  =  0: 


In  =  r,  where  Ryx{r)  =  0 


(4.34) 


Figure  4.13  contains  a  plot  of  In  for  all  £  after  t  =  5,  when  the  response  to  the  flow 
initial  condition  has  minimized.  Although  somewhat  noisy,  In  increases  in  time  and  is 
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Figure  4.13:  In  evolution  for  several  Note  the  horizontal  axis  begins  at  t  =  5,  when  the 
response  to  the  initial  condition  has  diminished. 

proportional  to  t°  5.  This  value  is  in  agreement  with  the  integral  length  scale  results  of 
Fraud  et  al.  [2005]  and  the  experimental  dipole  results  of  Voropayev  and  Afanasyev  [1992]. 


4.5.2  Energetics 
Kinetic  Energy 

A  general  idea  of  the  flow  energetics  can  be  seen  by  examination  of  the  x  direction  spectrum 
of  horizontal  energy  Efr(kx ),  defined  as: 

E/y ( kx )  —  ”  (u(kx)u(kx)  4*  v(kx}v(kx)  )  . 

Here  (•)  signifies  x-direction  Fourier  space,  (•)*  denotes  complex  conjugate,  and  kx  denotes 
x-direction  wavenumber.  The  left  plot  of  Figure  4.14(a)  contains  a  plot  of  En{kx)  for 
several  different  Rer  at  t  =  20.  The  initial  peak  at  kx  ~  tt/2  corresponds  to  the  separation 
distance  between  vortices.  The  large  scale  spectra  is  identical  for  all  Rer,  which  is  typical 
of  nearly  inviscid  large  scale  motions.  In  contrast,  the  small-scale  motions  (large  kx )  are 
shown  to  increase  monotonically  with  Rer,  consistent  with  the  ability  of  the  flow  to  generate 
small-scale  instabilities  with  increasing  Rer. 

A  time  evolution  of  En(kx)  for  Fr  =  2.75,  Rer  =  19200  is  shown  in  the  right  plot  of  Figure 
4.14(b).  As  the  flow  evolves,  small  scale  motions  are  generated  causing  the  magnitude  of 
En{kx)  at  high  wavenumbers  to  increase  between  t  =  0  and  t  =  10.  This  is  consistent  with 
the  increase  in  {S2)  ^  shown  in  Figure  4.11.  In  addition,  once  the  flow  has  become  turbulent 

(after  t  =  10),  it  displays  a  near  kx^3  spectrum.  This  kx  ^  dependence  is  seen  for  all  £ 
(Figure  4.15),  and  is  in  agreement  with  the  energy  spectra  results  of  Lindborg  [2006],  who 
also  found  a  /c~5/3  slope  for  horizontal  spectra  of  horizontal  energy.  After  t  =  10  the  energy 
at  all  wavenumbers  decreases  as  the  flow  decays. 
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Figure  4.14:  (a)  Ejj  x  spectra  for  several  Rer  at  t  —  10;  (b)  Time  evolution  of  Eh  x  spectra 
for  Fr  =  2.75,  Rer  =  19200.  Note  £  —  0.01  for  both  plots.  The  bold  dashed-dotted  line 
represents  kx 


Figure  4.15:  Eh  x  spectrum  for  all  t  =  10.  The  bold  dashed  line  represents  k  5/3.  The 
bold  dashed-dotted  line  represents  k~ 2 . 
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Figure  4.16:  Time  evolution  of  (a)  (Eh)  and  (b)  (Ey)  for  several  Rer,  Fr  =  2.75,  =  0.01. 


Figures  4.16  (a)  and  (b)  contain  plots  of  the  time  evolution  of  (Eh)  and  (Ey)  for  several 
Rer  and  £  =  0.01.  (Eh)  is  shown  to  be  relatively  independent  of  Rer.  In  eontrast,  the 
maximum  (Ey)  is  directly  related  to  Rer,  a  likely  result  of  the  increased  ability  for  vertical 
motion  to  overcome  the  stratification  as  Rer  is  increased.  Since  the  initial  vertical  velocity 
(ie)  and  fluctuating  density  (p)  are  set  to  zero,  the  flow  is  not  in  cyclostrophic  balance.  As 
a  result,  a  large  internal  wave  immediately  forms,  which  results  in  an  oscillation  in  both 
(Eh)  and  (Ey). 

The  time  evolution  of  (£*//)  and  (Ey)  for  each  £  is  shown  in  Figures  4.17(a),  (b).  As 
£  is  increased  from  0.01  to  2,  the  trend  is  for  (Eh)  to  persist  longer  in  time  and  the 
magnitude  of  (Ey)  to  decrease,  compared  with  the  unstratified  flow.  These  trends  are 
expected,  since  more  work  is  required  to  overcome  the  stronger  vertical  stratification  in 
the  vortex  “core.”  However,  these  trends  are  reversed  when  £  >  2,  and  the  flow  begins  to 
demonstrate  traits  similar  to  non-stratified  flow.  This  is  because  the  characteristic  height  of 
the  wake  is  much  smaller  than  the  density  stratification  height,  and  a  significant  portion  of 
the  wake  flow  occurs  in  regions  that  have  very  weak  density  stratification.  It  is  interesting 
that  the  transition  occurs  when  £  >  2,  rather  than  £  >  1,  suggesting  that  it  is  the  energy 
profile,  rather  than  the  velocity  profile,  that  determines  when  the  effect  of  stratification  is 
minimized.  Thus,  for  the  simulations  with  £  <  2,  the  stratification  acts  as  a  boundary  to 
confine  ( Ey)H .  In  eontrast,  when  (  =  4  a  significant  amount  of  (Ey)H  forms  in  regions 
outside  the  density  stratification  layer  where  there  is  no  buoyancy  force  and  overturning 
can  readily  occur.  This  can  be  seen  in  Figure  4.18,  where  (Ey) H  is  plotted  vs.  z  for  each 
£  and  no  stratification. 

T//  and  Ty  are  the  advection  (or  transport)  terms  that  describe  how  Ej)  and  Ey  arc 
adveeted  throughout  the  flow.  By  definition,  a  simulated  flow  with  periodic  boundary 
conditions  will  have  zero  advection  on  average,  since  any  transport  out  of  one  boundary 
will  be  transported  into  the  opposite  boundary.  Hence  (T//)  and  (Ty)  are  zero,  and  are  not 
shown  here. 

The  time  evolution  of  (T//) H  for  £  =  64)1  and  £  =  4  are  shown  in  Figures  4.19(a), 


4.5.  SIMULATION  RESULTS 


59 


Figure  4.17:  Time  evolution  of  (a)  {£//)  and  (b)  (Ey)  for  all  £. 


(b).  For  £  =  0.01,  early  in  the  simulation  when  turbulence  is  forming  (at  t  =  5)  there  is 
significant  advection  of  Ejj  into  the  center  planes  from  the  adjoining  ones  (recall  that  it  is 
negative  ( Tjj)u  in  (4.12)).  Advection  of  energy  into  the  center  planes  is  counter-intuitive, 
as  one  would  expect  energy  to  be  advected  from  the  maximum  energy  planes  to  minimum 
energy  planes.  This  can  be  explained  by  noting  that  the  planes  of  positive  (T//)/7  (z  ±0.5) 
are  also  planes  of  (S2)  H-  Since  planes  of  maximum  (S2)H  will  generate  a  large  amount  of 
vertical  motion  (through  shear  instabilities),  energy  will  be  advected  out  of  the  planes  of 
high  shear  into  the  lower  shear  center  plane.  (Tjj)  }J  most  diminishes  by  /  =  10. 

The  time  evolution  of  (Ty)ff  for  £  —  0.01  and  £  =  4  are  shown  in  Figures  4.20(a), 
(b).  Similar  to  (T//)/7,  there  is  a  large  amount  of  (Ty) H  early  in  the  flow.  However,  Ey 
is  advected  out  of  the  center  planes  into  the  adjoining  planes,  in  contrast  to  the  inward 
advection  of  E^.  The  reason  for  this  early  transport  out  of  the  center  of  Ey  is  unclear. 
Perhaps  since  vertical  overturning  is  originating  at  the  planes  of  maximum  shear,  eddies 
will  cause  a  net  advection  of  (Ey) n  out  of  the  center  planes. 

Pjl  and  Py  are  pressure  work  terms  that  describe  how  pressure  gradients  affect  the 
flow  energetics.  When  the  dot  product  of  velocity  and  the  pressure  gradient  taken,  the 
result  is  a  pressure  work  term  that  affects  the  evolution  of  energy.  For  example,  a  pressure 
gradient  in  the  direction  of  the  flow  increases  energy,  while  flow  opposite  a  pressure  gradient 
expends  energy.  Volume  averaged  (Pji)  and  (Py)  for  several  Rer  and  £  =  0.01  are  shown 
in  Figures  4.21(a),  (b).  Initially  large  pressure  work  is  induced  as  a  result  of  internal  wave 
generation  from  the  initial  cyclostrophic  imbalance.  This  early  pressure  work  results  in  an 
initial  loss  of  Efj  and  gain  in  Ey.  Since  (Py)  is  initially  negative,  when  substituted  into 
(4.13),  it  will  cause  an  increase  in  Ey.  The  opposite  is  true  for  (Pji)  and  (Eu).  Note  that 
for  a  closed  system  the  pressure  can  do  no  net  work,  so  that  the  total  pressure  work  term 
(P)  =  (Pii)  +  (Py)  is  zero.  It  is  interesting  that  the  magnitudes  of  (Ph)  and  (Py)  are 
nearly  identical  for  each  Rer,  suggesting  internal  wave  motion  is  controlled  by  large  scale 
dynamics. 

Figure  4.22(a)  contains  a  plot  of  (Pa)  vs.  time  for  each  £.  Figure  4.22(b)  contains  the 
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Figure  4.18:  (£y)//  f°r  several  £  and  no  stratification.  Note  the  difference  in  horizontal 
scale  for  plots  (e),  (f).  dp(z)/dz  is  shown  as  bold  dashed  line,  and  is  scaled  differently  in 
each  figure  for  comparison  to  (Ey)H.  Bold  horizontal  lines  mark  where  dp{z)/dz  =  —0.01. 
Since  the  entire  vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold  horizontal 
lines  are  shown. 
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Figure  4.19:  Evolution  of  (T//)  fI  for  (a)  £  =  0.01;  (b)  £  =  4.  Note  other  stratified  simulations 
are  qualitatively  similar 


Figure  4.20:  Evolution  of  (Ty) I{  for  (a)  £  =  0.01;  (b)  £  =  4.  Note  other  stratified  simulations 
are  qualitatively  similar 
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Figure  4.21:  (a)  (P//);  (b)  {Pv)  for  several  Rer  and  £  =  0.01 


same  data  as  Figure  4.22(a),  but  with  the  horizontal  axis  zoomed  in  to  t==  0-5.  Since  pressure 
work  must  sum  to  zero  for  a  periodic  simulation,  (P\  )  will  be  the  opposite  of  (P//)  (as  seen 
in  Figure  4.21  above)  and  is  not  shown  here.  Early  in  time,  there  is  a  significant  amount  of 
pressure  work  created  as  internal  waves  are  generated  in  response  to  the  initial  cyclostropliic 
imbalance.  It  is  interesting  to  see  that  the  magnitude  of  (P//)  for  £  =  4  is  nearly  half  the 
magnitude  of  (P//)  for  the  other  £.  The  reduction  in  magnitude  for  £  —  4  can  be  explained 
by  noting  that  internal  waves  cannot  exist  when  there  is  no  density  stratification,  and  when 
(  =  4  a  significant  portion  of  the  simulation  domain  has  minimal  density  stratification 
(Figure  4.7).  Thus,  internal  waves  (and  hence  and  pressure  work)  do  not  occur  in  much  of 
the  simulation  area  for  £  =  4,  and  a  drop  in  the  volume  averaged  (P//)  is  seen.  This  can  be 
seen  in  a  plot  of  the  horizontally  averaged  (P//)//,  shown  in  Figure  4.23. 

The  kinetic  energy  dissipation  rates  £//  and  ey  are  the  amount  of  E\\  and  Ey  that  are 
dissipated  into  internal  energy  via  viscous  effects.  The  time  evolution  of  (eu)  and  (ey) 
with  £  =  0.01  and  several  different  Rer  are  shown  in  Figures  4.24(a),  (b).  Since  the  initial 
conditions  are  the  same  for  each  simulation,  and  only  Reynolds  number  is  changed,  the 
initial  dissipation  rates  (eh)  and  (ey)  are  inversely  proportional  to  Rer.  This  relationship 
changes  as  the  flow  develops;  (eu)  increases  with  Rer  as  the  flow  has  more  ability  to  generate 
small-scale  instabilities,  consistent  with  the  increase  in  small-scale  spectra  shown  in  Figure 
4.14(a).  It  is  interesting  to  note  that  while  (£//)  and  (ey)  increase  faster  with  Rer,  their  final 
magnitudes  are  nearly  the  same  for  all  Rer.  This  suggests  that  a  “high  Reynolds  number 
limit”  has  been  reached,  where  the  flow  cannot  further  dissipate  kinetic  energy. 

Copmarisons  of  (eh)  and  (ey)  for  each  £  simulation  are  found  in  Figures  4.25(a),  (b). 
As  the  flows  evolve,  small  scale  turbulence  forms  and  both  (eu)  and  (ey)  increase  in  time 
before  decaying.  There  is  a  clear  transition  in  the  behavior  of  (£//)  and  (ey)  between  £  =  2 
and  £  =  4.  The  reason  for  this  transition  can  be  explained  by  the  behavior  of  (S2)  H,  the 
terms  of  which  are  components  of  (eh)  and  (ey).  When  £  =  4,  a  large  amount  of  overturning 
occurs  outside  the  density  stratification  layer,  which  causes  a  significant  amount  of  shear  and 
kinetic  energy  dissipation  rate  to  form  early  in  the  simulation  (Figures  4.26,  4.27.  Figure 
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Figure  4.22:  (a)  Time  evolution  of  (Ph)  for  several  Fr  =  2.75.  Re,  =  19200.  (b)  Same  as 
(a),  but  zoomed  horizontal  axis. 


Figure  4.23:  Evolution  of  {Ph)h  for  (a)  £  =  0.01;  (b)  £  =  4. 
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Figure  4.24:  Evolution  of  (a)  (£//),  (b)  (ey)  for  all  £  =  0.01  simulations. 


4.28  contains  plots  of  (ep)H  and  ( sy)fI  for  £  =  4  with  zoomed  vertical  axis). 

0  is  the  quantity  that  remains  when  separating  the  last  term  in  the  total  kinetic  energy 
equation  (4.6)  into  horizontal  and  vertical  components.  It  is  a  term  that  couples  horizontal 
and  vertical  kinetic  energy  through  a  process  of  viscous  and  rotational  strain.  The  time 
evolution  of  ((f))  for  all  £  =  0.01  simulations  is  shown  in  Figure  4.29(a).  Similar  to  (ep),  ((f)) 
increases  with  turbulence  formation  lip  to  approx,  t  =  8.  However,  in  contrast  with  (ep) 
and  ( ey ),  there  is  a  clear  inverse  relation  between  ((f))  and  Re  throughout  the  simulation. 
This  is  likely  due  to  the  decrease  in  v  with  Rcr,  causing  a  decrease  in  viscous  coupling 
(and  hence  0)  between  (Ep)  and  (Ey).  Note  the  magnitude  of  (f)  is  much  smaller  than  e, 
indicating  a  very  small  role  in  the  transport  of  energy. 

Figure  4.29(b)  contains  a  plot  of  ((f))  for  each  £,  where  the  overall  trend  is  for  ((f))  to 
increase  with  £  up  to  £  =  4.  The  increase  in  magnitude  of  (0)  with  £  can  be  explained 
by  noting  that  as  £  is  increased  the  local  stratification  is  increased.  The  increase  in  local 
stratification  causes  “stronger”  decoupling  of  horizontal  layers,  and  an  increase  in  shear 
that  results  in  stronger  coupling  between  vertical  and  horizontal  motions.  When  £  =  4  ((f)) 
reaches  a  maximum  much  earlier  than  the  other  simulations.  This  behavior  is  similiar  to 
that  seen  for  (ep)  and  (ey)  (Figure  4.25),  where  strong  horizontal  decoupling  causes  a  large 
amount  of  (0)  to  form  early  in  the  simulation. 

It  is  interesting  that  in  Figure  4.29  (0)  is  nearly  zero  for  the  non- stratified  case.  This 
seems  counter-intuitive,  since  when  there  is  no  stratification  there  is  no  buoyancy  force  and 
overturning  can  readily  occur.  However,  when  there  is  no  stratification  there  is  very  little 
horizontal  layer  decoupling  compared  to  stratified  flow.  With  no  buoyancy  force,  rotational 
shear  is  minimized,  and  0  will  be  much  smaller  compared  to  when  the  flow  is  stratified.  This 
is  seen  in  plots  of  the  planar  averaged  (0) u  shown  Figure  4.30,  where  ((f)) H  is  barely  visible 
for  the  No  Strat  simulation.  Also,  when  £  =  4  ((f)) H  is  near  zero  outside  the  stratification 
layer,  whereas  significant  (e) H  occurs  when  the  stratification  is  minimal  (Figures  4.26  and 
4.27). 

The  buoyancy  flux  term  B  is  a  coupling  term  between  Ey  and  Ep.  As  mass  is  advected 
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Figure  4.25:  Evolution  of  (a)  (£//),  (b)  (ey)  for  each  £  simulation. 


vertically  it  must  do  work  against  gravity,  and  energy  is  transferred  from  Ey  to  Ep.  F  igure 
4.31(a)  contains  the  time  evolution  of  (B)  for  several  Rer,  £  =  0.01.  Interpretation  of  (B) 
on  a  time  basis  is  difficult  because  Ey  converted  to  Ep  at  one'  time  can  be  released  back  to 
Ey  at  another.  However,  from  Figure  4.31(a)  one  can  conclude  that  there  is  a  large  initial 
buoyancy  flux  (B)  as  the  flow  adjusts  to  the  initial  cyelostrophic  imbalance,  supported 
by  the  fact  that  the  oscillation  in  (B)  is  exactly  180°  out  of  phase  with  ( Py ).  After  this 
initial  adjustment,  the  remaining  oscillations  in  ( B )  are  likely  caused  by  the  propagation  of 
internal  gravity  waves. 

Potential  Energy 

Figure  4.32(a)  contains  a  plot  of  the  volume  averaged  available  potential  energy  (Ep)  for 
several  Rer,  Fr  =  2.75,  £  =  0.01,  where  a  rapid  initial  increase  in  (Ep)  is  caused  by  the 
formation  of  an  internal  wave  from  initial  cyelostrophic  adjustment.  Following  this  initial 
rise,  (Ep)  decays  to  an  undisturbed  state.  The  oscillations  in  (Ep)  after  t  ~  5  are  transfer 
between  Ep  and  Ey  via  buoyancy  flux  (B).  This  is  supported  by  the  the  fact  that  the 
oscillations  in  (Ej>)  and  (B)  (Figure  4.31(a))  are  in  phase  with  each  other.  The  increase  in 
(Ep)  with  Rer  is  due  to  the  flow’s  increased  ability  to  overcome  stratification. 

The  time  evolution  of  (Ep)  for  each  £  is  show  in  in  Figure  4.32(b).  After  the  initial  flow 
response  to  the  initial  contitions  by  t  ss  5,  the  general  trend  is  for  (Ep)  to  be  relatively 
independent  of  £  until  £  =  2,  after  which  there  is  a  drop  in  (Ep)  with  £.  This  drop  in  (Ep) 
can  be  explained  by  noting  that,  by  definition,  Ep  is  zero  outside  the  stratification  layer 
(since  in  the  limit  where  dp(z)/dz  =  0,  N  =  0,  and  Fr  =  oo).  Recalling  that  wp  is  a  source 
term  of  Ep  (equation  4.26),  when  £  =  2,  the  vertical  span  of  Ep  has  narrowed  beyond  the 
vertical  span  of  w.  (Ep)  will  thus  incorporate  areas  where  Ep  is  zero,  even  though  vertical 
motion  occurs.  This  can  be  seen  in  Figure  4.33,  where  the  peak  planar  averaged  (Ep)  H  is 
nearly  the  same  for  all  £,  even  though  there  a  sharp  increase  in  the  volume  averaged  (Ep). 

Tp  is  a  term  that  describes  the  advection  of  potential  energy,  and  is  similar  to  T//  and 
Ty.  Since  the  boundary  conditions  are  periodic,  the  volume  averaged  (Tp)  will  be  zero,  and 
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Figure  4.26:  Evolution  of  (sh)h  for  all  The  dashed  line  represents  the  density  stratifi¬ 
cation  dp(z)/dz.  The  bold  solid  lines  identify  where  dp(z)jdz  =  —0.01.  Since  the  entire 
vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold  horizontal  lines  are  shown. 
When  £  =  4,  much  of  the  (sh)h  is  formed  earlier  in  time  outside  the  stratification  layer 
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Figure  4.27:  Evolution  of  {sy)  n  for  all  £.  The  dashed  line  represents  the  density  stratifi¬ 
cation  dp(z)/dz.  The  bold  solid  lines  identify  where  dp(z)/dz  =  —0.01.  Since  the  entire 
vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold  horizontal  lines  are  shown. 
When  £  =  4,  much  of  the  (sy)lf  is  formed  earlier  in  time  outside  the  stratification  layer 
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Figure  4.28:  (e//)//  (left)  and  (ev) //  (right)  for  £  =  4.  The  plots  are  the  same  as  in  Figures 
4.27(e)  and  4.27(e),  but  with  zoomed  vertical  axis. 


Figure  4.29:  Evolution  of  (<f>)  for  (a)  several  Rer,  £  =  0.01,  (b)  Several  £,  Fr  =  2.75, 
Rer  =  19200.  The  lines  in  (b)  correspond  to  those  shown  in  Figure  4.25 
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Figure  4.30:  (</>)/./  fit,  for  all  £  at  t  =5  (left)  and  t  =10  (right).  Note  that  when  the 
stratification  is  small  for  £  =  4  and  the  No  Strat  simulation  (</>)w  is  small. 


Figure  4.31:  Evolution  of  ( B )  for  (a)  several  Rer.  £  =  0.01,  (b)  Several  £,  Fr  =  2.75, 
Rer  =  19200. 
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Figure  4.32:  Time  evolution  of  ( Ep )  for  (a)  several  Rer,  Fr  =  2.75,  f  =  0.01,  (b)  Several  f, 
Fr  =  2.75,  Rer  =  19200. 

is  not  shown  here. 

Planar  averaged  (Tp)H  for  F2.75R192  is  shown  in  Figure  4.34.  (Tp)n  is  qualitatively 
similar  for  all  simulations.  Early  in  the  simulation  at  t  ^  5,  Ep  is  seen  to  be  advected  into 
the  center  planes  from  the  adjoining  planes,  where  the  initial  high  vertical  shear  is  formed. 
Advection  then  diminishes,  as  the  majority  of  transport  of  Ep  occurs  due  to  D  and  x- 

X  is  the  potential  energy  dissipation  rate,  which  can  be  described  as  the  irreversible 
conversion  from  available  to  background  potential  energy,  similar  to  e  (see  §4.2.2  for  further 
discussion  on  background  potential  energy).  The  time  evolution  of  the  volume  averaged  (x) 
for  several  Rer  and  £  =  0.01  is  shown  in  Figure  4.35(a).  Since  the  flow  is  initialized  with  no 
density  fluctuation,  (x)  is  also  initially  zero.  Early  in  the  simulation  (t  <  7)  (x)  is  inversely 
proportional  to  Rer,  a  result  of  the  fact  that  the  initial  condition  for  each  simulation  is  the 
same,  and  only  the  mass  diffusivity  V  is  changed  (recall  that  Sc  =  1  for  all  simulations). 
Then,  as  the  flow  evolves  further  and  turbulent  small-scale  motions  form,  (x)  demonstrates 
an  increase  with  Rer,  similar  to  (e). 

The  time  evolution  of  (x)  for  each  £  is  shown  in  Figure  4.35(b),  where  the  magnitude  of 
(x)  increases  in  time  as  small  scale  turbulence  forms.  The  time  to  reach  peak  (x)  increases 
as  £  is  increased  from  0.01  to  2.  When  £  =  4,  in  contrast,  the  peak  (x)  occurs  earlier, 
indicating  turbulent  motions  occur  earlier  in  time.  This  behavior  is  in  agreement  with  that 
seen  for  the  kinetic  energy  dissipation  rates  (Figure  4.25).  However,  unlike  (e//)  and  (sy), 
the  magnitude  of  (x)  decreases  when  £  >  2.  This  is  because  x  is  zero  when  there  is  no 
stratification  for  the  same  reason  Ep  is  zero.  Thus,  when  £  >  2,  a  significant  portion  of  the 
volume  has  x  =  0,  (Figure  4.36),  which  causes  (x)  to  be  small. 

4.5.3  Mixing  and  Mixing  Efficiency 

The  behavior  of  (e)  and  (x)  with  respect  to  £  has  implication  in  modeling  the  Earth’s 
energy  budget.  Current  Earth  general  circulation  models  (gem’s),  assuming  linear  density 
stratification,  show  an  imbalance  in  vertical  heat  fluxes  by  as  much  as  20%  [Howard  et  al., 
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Figure  4.33:  Evolution  of  ( Ep)fJ  for  each  dp(z)/dz  is  shown  as  bold  dashed  line,  and  is 
scaled  differently  in  each  figure  for  comparison.  Bold  horizontal  lines  mark  where  dp(z)/dz  = 
—0.01.  Since  the  entire  vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold 
horizontal  lines  are  shown. 
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Figure  4.34:  Evolution  of  {Tp)u  for  simulation  F2.75R192 


Figure  4.35:  Evolution  of  (\)  for  (a)  several  Rer,  £  =  0.01,  (b)  Several  £,  Fr  =  2.75, 
Rer  =  19200. 
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Figure  4.36:  Evolution  of  (\)  for  each  dp(z)/dz  is  scaled  different  for  each  figure,  and  is 
shown  for  comparison  to  (x)h*  Bold  horizontal  lines  mark  where  dp(z)/dz  =  —0.01.  Since 
the  entire  vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold  horizontal  lines 
are  shown. 
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2004,  Marmorino  and  Caldwell,  197G,  Robertson  et  ah,  1995].  This  imbalance  may  be  due 
to  the  mixing  results  obtained  with  a  linear  density  stratification  throughout  the  Earth’s 
atmosphere  and  ocean.  From  Figures  4.25  and  4.35(b),  one  can  see  that  for  the  same  initial 
flow  conditions,  peak  (e)  increases  by  25%  when  a  sharp  density  gradient  exists  (£  =  4), 
while  (x)  decreases  by  50%.  This  results  suggests  that  in  areas  of  sharp  density  gradients, 
assuming  a  uniform  density  gradient  leads  to  under-predicting  mixing  of  kinetic  energy  e, 
and  over-predicting  x* 

As  discussed  in  §1.2.2,  a  quantity  of  interest  is  the  mixing  efficiency  T,  which  is  the  ratio 
of  the  available  potential  energy  lost  by  molecular  diffusion  to  kinetic  energy  lost  by  viscous 
dissipation  [e.g.,  Smyth  et  ah,  2001,  Winters  et  ah,  1995].  Here  T  will  be  defined  as  in  (1.4) 
(and  dropping  the  subscript  a): 

r  =  |.  (4.35) 

F  is  of  interest  because  it  is  difficult  to  measure  both  e  and  x  simultaneously  in  field 
experiments  due  to  the  time  scales  required  to  measure  each  quantity  [Gargett  and  Mourn, 
1995].  Thus,  it  would  be  beneficial  to  measure  one  quantity  (usually  e)  and  relate  it  to 
the  other  [Osborn,  1980].  Figure  4.37(a)  contains  the  volume  averaged  mixing  efficiency 
(F)  =  (x)  /  (?}•  Initially,  F  =  0  (since  x  =  0)*  after  which  there  is  a  rapid  increase  in  (r) 
to  near  0.9  as  the  density  field  responds  to  the  initial  cyclostrophic  imbalance.  As  the  flow 
evolves  and  small-scale  instabilities  occur,  e  increases  and  causes  (F)  to  drop  to  ^  0.45,  in 
agreement  with  results  from  other  numerical  simulations  of  uniformly  stratified  flows  [e.g., 
Riley  and  de  Bruyn  Kops,  2003  Smyth  et  ah,  2001.  Staquet,  2000].  (F)  is  also  shown  to  be 
relatively  independent  of  Rer,  which  is  consistent  with  the  aforementioned  studies. 

The  evolution  of  (Y)  for  each  £  is  shown  in  Figure  4.37(b),  where  the  general  trend  is  for 
a  large  initial  (Y)  to  form  up  to  t  ^  7,  after  which  (F)  settles  to  a  constant  value  related  to 
the  value  of  £.  The  early  rise  in  (r)  is  due  to  the  flow’s  initial  ability  to  dissipate  available 
potential  energy  to  background  potential  energy  (x)  faster  than  kinetic  energy  is  dissipated 
to  internal  energy  (f).  As  the  flow  progresses  and  small-scale  instabilities  occur,  e  increases 
and  causes  (F)  to  drop  before  settling  to  some  constant  value.  These  results  are  consistent 
with  the  aforementioned  studies  of  uniformly  stratified  flows. 

Since  the  vertical  span  of  x  decreases  as  £  increases  (Figure  4. 36),  the  value  of  (F) 
will  also  decrease  due  to  the  inclusion  of  areas  of  zero  mixing  in  the  volume  average.  It 
would  be  of  interest  to  know  what  the  mixing  efficiency  is  inside  the  density  stratification 
layer.  In  this  case,  the  density  stratification  layer  is  defined  where  \dp(z)/dz\  >  0.01.  The 
value  of  0.01  is  chosen  because  it  is  1%  of  the  stratification  value  when  £  =  0.01.  Figure 
4.38(b)  contains  a  plot  of  the  time  evolution  of  {Y)STRAT  =  ( x)strat  /  (£) strata  w^iere 
the  subscript  ST  RAT  denotes  average  over  the  density  stratification  range.  Large  initial 
values  of  ( Y )strat  can  be  seen  f°r  the  same  reasons  as  (r).  It  is  interesting  to  note  that 
(F) st RAT  set^es  to  near  the  same  value  for  all  £. 

To  get  a  better  idea  of  how  F  is  behaving  throughout  this  flow,  the  planar  averaged 
(r) H  =  (x) n  /  (e) h  for  each  £  is  shown  in  Figure  4.39,  where  the  value  of  (T) H  is  relatively 
independent  of  £  in  the  center  of  the  domain.  The  peak  at  2  =  0  is  due  to  the  centerplane 
dip  in  (e)H  seen  in  Figures  4.26,  4.27.  The  trend  for  (Y)  H  to  increase  when  \z\  >  0.5  is  due 
to  the  decreasing  value  of  e  in  this  range.  Outside  \z\  >  0.5  e  is  very  small  and  results  in  a 
noisy  {Y)H,  particularly  for  f  =  0.5.  In  addition,  since  x  is  zero  outside  the  stratification 
layer,  Y  is  also  zero.  This  suggests  that  in  an  area  with  non-uniform  density  stratification, 
mixing  can  only  occur  in  the  region  of  high  density  gradient.  Such  a  result  is  in  agreement 
with  the  findings  of  Schmitt  [2003],  St.  Laurent  and  Schmitt  [1999],  where  increased  mixing 
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Figure  4.37:  (F)  for  (a)  several  Rer,  £  =  0.01;  (b)  all  Fr  =  2.75,  Rer  =  19200. 


Figure  4.38:  Evolution  of  (L)STRAT  for  each 
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is  found  in  the  sharp  density  gradients  of  therinohaline  staircases,  where  salt  fingering  is 
prone  to  occur. 


4.5.4  Buoyancy  Reynolds  Number 


In  stratified  Hows,  the  importance  of  turbulence  is  often  parameterized  in  terms  of  the 
buoyancy  Reynolds  number  Reb  [e.g.,  Gibson.  1980,  Gregg,  1987,  Imberger  and  Boashash, 
1986,  Smyth  and  Mourn,  2000a].  Reb  is  usually  defined  as  the  ratio  of  the  Ozmidov  scale  L0 
to  the  Kolmogorov  scale  (see  §2.1.9  for  further  discussion  on  Reb).  Recalling  equation 
2.57,  Reb  is  defined  as: 


Reb  = 


(4.36) 


From  this  definition,  Reb  =  21.5  for  one  decade  of  turbulent  length  scales  to  form.  Analysis 
by  Gibson  [1980]  resulted  in  Reb  ~  30  for  “active”  turbulence  to  form.  Figure  4.40(a)  is 
a  plot  of  the  time  evolution  of  the  volume  averaged  buoyancy  Reynolds  number  (Reb)  = 

(e)  j(v N  )  for  several  Rer,  where  the  maximum  (Reb)  of  14  might  suggest  that  the  flow  does 
not  have  an  adequate  range  of  length  scales  for  turbulence  to  form.  This  is  misleading,  as 
the  volume  averaged  data  includes  areas  outside  the  wake  “core”  where  there  is  little  energy 
(and  hence  ?),  particularly  at  the  outer  edges  of  the  flow  where  \z\  >  2.  Consequently,  the 

planar  averaged  (Reb) H  —  (e) u  jv  (n ^  can  be  examined  to  get  a  better  idea  of  how  Reb  is 

acting  locally.  Figure  4.40(b)  contains  a  plot  of  the  time  evolution  of  (Reb)//  for  simulation 
F2.75R192,  where  the  magnitude  of  (Reb)//  approaches  40  in  the  center  of  the  vertical 
domain  by  t  =  10.  Thus,  according  to  the  criterion  specified  by  Gibson,  the  simulated  flow 
is  in  fact  actively  generating  turbulent  motions.  It  is  also  interesting  to  note  the  dip  in 
(Reb)//  at  z  =  0,  which  occurs  at  the  same  location  as  the  dip  in  {S2) ^  (Figure  4.11), 
suggesting  a  strong  role  of  S2  in  the  formation  of  turbulent  motions. 

Figure  4.41(a),  contains  a  plot  of  the  volume  averaged  (Reb)  for  each  £.  Similar  to  Figure 
4.40(a),  the  maximum  value  of  (Reb)  is  less  than  the  21.5  needed  for  active  turbulence  to 
form.  Also,  at  t  =  10,  the  maximum  (Reb)//  is  between  22  and  40  (Figure  4.41(b)), 
indicating  active  turbulence  formation.  The  rapid  rise  in  (Rob)  for  £  =  4  is  attributed  to 
the  behavior  of  (e),  which  increases  rapidly  due  to  the  shear  caused  in  the  region  where 
overturning  flow  intersects  with  strongly  stratified  flow. 


4.5.5  Vertical  shear  vs.  kinetic  energy  dissipation  rate 

In  field  studies,  e  is  usually  inferred  from  measurements  of  one  or  two  components  of  the 
strain  rate  tensor  etj.  In  strong,  stable  density  stratification,  it  is  often  assumed  that  S2 
causes  most  of  £,  leading  to  the  following  relationship: 

uS2  k,  e.  (4.37) 

While  the  ratio  uS2 /e  «  0.9  has  been  shown  to  hold  for  a  large  range  of  Reynolds  numbers 
[Fincham  et  al.,  1996,  Praud  et  ah.  2005],  in  §3.4.1  of  this  document  the  above  relation  was 
found  to  be  true  only  when  Reb  <  0(1).  Figure  4.42  contains  plots  of  the  volume  averaged 

quantity  v  ^52^)  /  (e)  vs.  (Reb)  for  (a)  several  Rer,  £  =  0.01,  and  (b)  several  £,  Fr  =  2.75, 
Rer  =  19200.  In  each  plot,  the  heavy  dashed  line  signifies  the  value  i>  (S 2)  /(e)  —  0.26, 
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Figure  4.39:  Evolution  of  (L)H  for  each  dp(z)/dz  is  shown  as  bold  dashed  line,  and  is 
scaled  differently  in  each  figure  for  comparison.  Bold  horizontal  lines  mark  where  dp(z)/dz  ~ 
—0.01.  Since  the  entire  vertical  domain  satisfies  this  criterion  for  plots  (a)-(b),  no  bold 
horizontal  lines  are  shown. 
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Figure  4.40:  (a)  (Reb)  for  several  Rer.  (b)  (Ih'b)//  for  F2.75R192. 


Figure  4.41:  (a)  (Reb)  for  each  Fr  =  2.75,  Rer  =  19200.  (b)  (Reb)w  for  all  £,  t= 10. 


V  {  s2  }  /  (  £  } 
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Figure  4.42:  v  (^S2^  /  (£)  vs. 


Re^  for  several  Rer. 


which  is  the  approximate  value  expected  if  the  flow  were  isotropic  [Hinze.  1975].  As  in  §3.4.  1, 
(Reij)  increases  to  above  1  as  the  flow  evolves  and  turbulence  forms.  Correspondingly, 
/>^S2^}/  (e)  decreases  to  below  0.5,  signifying  S 2  accounts  for  less  than  half  the  kinetic 
energy  dissipation  rate.  (Note  that  with  the  influence  of  gravity  on  stratified  flows,  isotropy 
is  not  expected,  hence  v  (S2)  /  (e)  should  not  reach  the  isotropic  limit).  Then,  as  the  flow 
decays,  (Re^)  decreases,  and  the  ratio  v  (S2)  /  (s)  increases. 

In  order  to  see  if  (4.37)  holds  locally,  the  quantity  0  (^S2^  /  (e)H  is  shown  in  Figure 

4.43.  In  comparing  with  Figure  4.41(b),  one  can  see  that  locally,  v  ^S2^}  /  {£}//  ~  0-8  only 

for  £  =  4,  where  (Re\>)H  <  1. 
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Figure  4.43:  v  (s2^  /  (e)jj  for  each  £  at  t  =  10. 


Chapter  5 


Summary  and  Suggested  Future 
Work 


5.1  Summary 

5.1.1  Uniform  Density  Stratification 

A  series  of  high- resolution  direct  numerical  simulations  are  used  to  investigate  turbulence  in 
stably  stratified  flows.  The  simulated  flows  are  dominated  by  vortical  modes  and  decay  in 
time  since  there  is  no  input  of  energy  from  the  mean  flow.  In  this  regime,  the  Froude  number 
decreases  and  the  Reynolds  number  increases  in  time  so  that  a  strongly  stratified  flow  with 
turbulence  eventually  develops  except  when  the  initial  Reynolds  number  is  extremely  low. 
The  predominant  cause  of  turbulence  in  the  simulations  is  Kelvin-Helmholtz  instabilities 
that  result  when  the  flow  organizes  itself  into  quasi-horizontal  vortices  that  are  weakly 
coupled  in  the  vertical  direction. 

This  observation  that  shear  instabilities  are  the  primary  trigger  for  turbulence  in  the 
simulations  supports  the  theoretical  derivation  of  the  F^Re/,  parameterization  developed 
in  Riley  and  de  Bruyn  Kops  [2003],  namely  that  1/F^Re/j  is  related  to  the  Richardson 
number,  so  when  F^Re/j  >  (9(1)  turbulence  can  be  expected.  To  test  this  hypothesis,  the 
horizontal  length  scale,  L/j,  on  which  F^  and  Re/>  arc  based,  is  defined  as  the  correlation 
length  of  the  horizontal  velocity.  When  this  definition  is  used  it  is  observed  that  L^,  F 
and  Re/j  all  evolve  in  time  consistently  with  theoretical  predictions  for  all  of  the  simulation 
cases.  Furthermore,  l/F^Rc/i  ~  Ri  over  a  two  decade  span  of  Richardson  numbers.  This 
result  encourages  the  thought  that  F^Re^  can  be  used  a  priori  to  estimate  if  a  laboratory 
or  simulated  flow  will  involve  considerable  turbulence  provided  that  the  correlation  length 
of  the  horizontal  velocity  can  be  estimated  from  the  initial  and  boundary  conditions  for  the 
flow. 

FjRc*  is  compared  with  the  buoyancy  Reynolds  number  for  all  the  simulation  cases, 
and  it  is  found  that  the  two  quantities  are  nearly  the  same.  While  Reb  is  traditionally 
defined  as  the  ratio  of  the  overturning  length  scale  to  the  viscous  length  scale,  the  fact 
that  F^Re/i  ~  Reb  suggests  that  a  shear-based  argument  might  be  used  to  better  explain 
why  Reb  has  proven  to  be  a  useful  parameterization.  Also,  since  F^Re/j  involves  the  num¬ 
ber  of  dimensionless  groups  predicted  by  dimensional  analysis,  it  encourages  considering 
turbulence  in  stratified  flows  as  occurring  when  the  Froude  and  Reynolds  numbers  are  in 
some  region  of  the  two-dimensional  Froude-Reynolds  number  space  rather  than  when  the 
buoyancy  Reynolds  number  is  above  some  transition  value. 
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Finally,  it  is  found  that  the  contribution  of  vertical  shear  to  the  total  dissipation  rate 
of  kinetic  energy  is  a  strong  function  {Reb).  The  simulation  results,  taken  in  conjunction 
with  the  laboratory  experiments  of  Fincham  et  ah  [1996]  and  Fraud  et  al.  [2005],  indicate 
that  the  relation  v  (S2)  /  (e)  ~  0.9  only  applies  when  (Reb)  is  order  one  or  less.  For  higher 
values  of  (Reb),  normal  strains  make  a  significant  contribution  to  (e)  and  (S2)  is  not  well 
correlated  with  (e).  As  a  result,  field  measurement  of  only  one  or  two  components  of  the 
strain  rate  in  the  ocean  can  lead  to  under-prediction  of  the  true  dissipation  rate  by  as  much 
as  a  factor  of  7. 


5.1.2  Nonuniform  Density  Stratification 

Simulations  are  conducted  for  Hows  initialized  with  a  von  Harman  vortex  street  and  no  mean 
velocity  or  shear,  similar  to  a  rnomentuinless  wake.  Each  simulation  was  subject  to  a  non- 
uniform  density  stratification  which  resembles  the  stratification  found  in  natural  settings 
such  as  a  therrnohaline  staircase.  A  method  is  derived  to  determine  available  potential 
energy  for  this  particular  flow,  followed  by  analysis  of  flow  energetics.  For  each  simulation, 
the  overall  change  in  density  with  height  is  the  same,  but  the  wake  height  relative  to  the 
density  layer  is  altered.  The  results  are  also  compared  to  a  simulated  flow  in  which  no 
density  stratification  is  present. 

The  simulated  density  stratified  flows  in  which  the  wake  height  is  less  than  or  equal  to 
twice  the  density  layer  height  (£  <  2)  are  shown  to  agree  with  the  current  understanding 
of  density  stratified  flows.  This  includes  growth  of  horizontal  length  scales  and  decrease 
in  vertical  length  scales  as  the  flow  evolves,  in  agreement  with  the  scaling  arguments  of 
Riley  et  al.  [1981]  and  the  horizontal  layer  decoupling  heuristic  by  Lilly  [1983].  However, 
when  the  wake  is  greater  that  twice  the  density  layer  height  (£  >  2),  the  importance  of 
the  density  stratification  is  diminished,  and  the  flow  begins  to  demonstrate  characteristics 
of  non-stratified  flows.  That  is,  Ep  dissipates  faster  and  Ey,  e,  and  x  demonstrate  rapid 
increases  in  magnitude.  The  transition  point  of  £  >  2  suggests  that  it  is  the  relation  of  the 
local  energy  profile  to  the  density  layer,  rather  than  velocity  profile,  that  will  determine  if 
the  flow  will  behave  in  a  stratified  manner. 

The  results  also  demonstrate  that  when  a  sharp,  localized  density  gradient  exists,  as¬ 
suming  a  uniform  density  gradient  will  under-report  (e)  up  to  25%,  while  over-reporting  (x) 
by  as  much  as  50%.  This  may  help  explain  errors  in  the  reporting  of  energy  budgets  from 
general  circulation  models  currently  in  use.  In  addition,  it  has  been  shown  that  available 
potential  energy,  Ep ,  and  potential  energy  dissipation  rate,  x>  are  confined  to  the  area  of 
density  stratification.  Thus,  when  £  >  2,  the  volume  averaged  (Ep)  drops,  while  locally 
(Ep) h  similar  in  magnitude  for  all  £. 

Finally,  the  mixing  efficiency  ( T )  is  shown  to  drop  significantly  when  £  >  2.  However, 
since  x  is  zero  when  there  is  no  stratification,  the  average  mixing  efficiency  inside  the 
stratification  region  Y STRAT  was  computed.  Y st RAT  was  found  to  be  between  0.3  and  0.5 
for  all  simulations,  which  is  in  agreement  with  previous  studies  of  uniformly  stratified  fluids. 
Also,  the  notion  that  Y  will  be  confined  to  the  area  of  density  stratification  by  definition 
lends  support  to  the  results  of  St.  Laurent  and  Schmitt  [1999],  in  which  large  values  of  Y 
are  reported  in  regions  where  salt  fingering  occurs. 
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5.2.1  Scales  of  Motion 

Turbulent  flows  subject  to  strong  stable  stratification  evolve  differently  from  flows  in  which 
the  stabilizing  effect  of  gravity  is  cither  absent  or  negligible.  Two  reasons  for  this  are  (1) 
potential  energy  plays  an  important  role  in  stratified  flows  and  (2)  stably  stratified  flows 
tend  to  be  quasi-two-dimensional  and  hence  lose  less  kinetic  energy  to  heat  through  viscous 
mechanisms  than  similar  mist  ratified  flows.  While  these  two  differences  between  stratified 
and  unstratified  turbulence  are  readily  apparent,  just  how  they  and  other  phenomena  affect 
the  flow  dynamics,  and  the  implications  with  respect  to  our  ability  to  predict  stably  stratified 
turbulence,  is  not  well  understood.  For  example,  questions  remain  as  to  what  scales  energy  is 
converted  from  one  form  to  another,  c.g.,  from  potential  to  kinetic,  kinetic  to  heat,  etc.  Also, 
there  is  a  question  of  how  the  scales  of  motion  are  affected  when  exposed  to  different  density 
stratifications.  The  above  questions  can  be  examined  using  direct  numerical  simulations, 
since  all  terms  in  the  momentum  and  density  equations  are  directly  computed,  and  all 
dynamically  relevant  scales  of  the  flow  are  adequately  resolved.  (This  is  not  to  say  that 
DNS  can  be  used  for  large  scale  simulations,  for  which  the  required  spatial  resolution  to 
perform  a  DNS  makes  such  simulations  impractical.  Rather,  DNS  is  a  powerful  tool  when 
used  for  the  appropriate  sized  case.)  It  is  also  of  interest  to  compare  how  flow  energetics 
differ  (if  at  all)  between  the  Taylor-Green  and  vortex  street  simulations. 

5.2.2  Turbulent  Patch  Identification  and  Tracking 

Turbulence  in  density  stratified  flows  tends  to  form  in  intermittent  patches,  and  requires 
a  large  number  of  data  points  to  adequately  describe  the  flow  characteristics  [Baker  and 
Gibson,  1987,  Gibson,  1981].  Also,  it  is  in  these  turbulent  patches  that  the  majority  of 
quantities  of  interest  occur,  including  mixing  and  dissipation  of  energy.  The  ability  to  iden¬ 
tify  and  track  turbulent  patches  could  lead  to  further  insight  into  fundamental  turbulent 
processes.  Attempts  at  identifying  turbulent  structures  have  been  made  by  using  wavelet 
analysis  [e.g.,  Dallman  et  ah,  1999,  Fargo,  1992,  Katul  and  Vidakovie,  1998],  which  is  sim¬ 
ilar  to  the  filtering  process  of  large  eddy  simulations  (LES).  The  difference  is  that  while 
LES  separates  the  fields  in  large  and  small  scales,  wavelet  analysis  separates  the  field  into 
Gaussian  and  non-Gaussian  components  [Goldstein  et  ah,  2000].  The  Gaussian  component 
represents  random  noise,  while  the  non-Gaussian  components  represents  coherent  struc¬ 
tures.  A  drawback  of  this  method  is  that  the  portion  of  the  flow  identified  a s  a  turbulent 
structure  is  dependent  on  the  type  of  wavelet  filter  used  (e.g..  Haar,  Donoho,  etc.)  and  the 
level  of  filtering  applied,  leaving  the  selection  of  coherent  structures  a  subjective  process. 
A  more  objective  method  would  remove  much  of  the  arbitrariness  involved  in  the  wavelet 
method,  and  perhaps  yield  clearer  turbulent  phenomena  insight. 
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Appendix  A 

Equation  of  state  for  seawater 


The  equation  of  state  for  the  density  of  seawater  (p)  is  given  in  terms  of  potential  tem¬ 
perature  (B0*?),  practical  salinity  (5,  in  practical  salinity  units),  and  pressure  ( p ,  in  bars) 
[Millero  and  Poisson,  1981]: 

p(S,  @,p)  =  p{S,  ©,0)/[l  —  pK(S,Q,p)} ,  (A.l) 

where  p(5,0,  0)  is  the  density  at  reference  pressure  of  zero  bars,  and  and  K  is  the  secant 
bulk  modulus: 

p(5.  0, 0)  =  999.842594  +  6.793952  x  10  ‘2©  -  9.09529  x  10 ~302 

+  1.001685  x  10  4©3  -  1.120083  x  10-6©4  +  6.536332  x  10-9©5 
+8.24493  x  10  lS  -  4.0899  x  10“35©  +  7.6438  x  10_5©25 
-8.2467  x  10  7©35  +  5.3875  x  10_9©4S  -  5.72466  x  10~3515 
+  1.0227  x  10  l©5‘  5  -  1.6546  x  lO'6©2^1  5  +  4.8314  x  10“  'S2 


I\(S,  e,p)  =  19652.21  +  148.4206©  -  2.327105©2  + 

1.360477  x  102©3  -  5.155288  x  105©4  +  3.239908p 
+  1.43713  x  103©p+  1.16092  x  104  ©2p  -  5.77905  x  107©3/~> 
+8.50935  x  10V  -  6.12293  x  10G©p2  +  54.67465 
-0.603459©5  +  1.09987  x  102©  25  -  6.167  x  105©35 
+7.944  x  102515  +  1.6483  x  102©5X  5  -  5.3009  x  104©25!  5 
+2.2838  x  103p5-  1.0981  x  105©/35  -  1.6708  x  106©2p5 
+  1.91075  x  lO4^1  5  -  9.9348  x  107p25  +  2.0816  x  108©p25 
+9.1697  x  1010©2p25; 
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